1 #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/ 2 3 /*@ 4 DMNetworkGetPlex - Gets the Plex DM associated with this network DM 5 6 Not collective 7 8 Input Parameters: 9 . dm - the dm object 10 11 Output Parameters: 12 . plexdm - the plex dm object 13 14 Level: Advanced 15 16 .seealso: DMNetworkCreate() 17 @*/ 18 PetscErrorCode DMNetworkGetPlex(DM dm,DM *plexdm) 19 { 20 DM_Network *network = (DM_Network*)dm->data; 21 22 PetscFunctionBegin; 23 *plexdm = network->plex; 24 PetscFunctionReturn(0); 25 } 26 27 /*@ 28 DMNetworkGetNumSubNetworks - Gets the the number of subnetworks 29 30 Not collective 31 32 Input Parameters: 33 . dm - the dm object 34 35 Output Parameters: 36 + nsubnet - local number of subnetworks 37 - Nsubnet - global number of subnetworks 38 39 Level: beginner 40 41 .seealso: DMNetworkCreate(), DMNetworkSetNumSubNetworks() 42 @*/ 43 PetscErrorCode DMNetworkGetNumSubNetworks(DM dm,PetscInt *nsubnet,PetscInt *Nsubnet) 44 { 45 DM_Network *network = (DM_Network*)dm->data; 46 47 PetscFunctionBegin; 48 if (nsubnet) *nsubnet = network->nsubnet; 49 if (Nsubnet) *Nsubnet = network->Nsubnet; 50 PetscFunctionReturn(0); 51 } 52 53 /*@ 54 DMNetworkSetNumSubNetworks - Sets the number of subnetworks 55 56 Collective on dm 57 58 Input Parameters: 59 + dm - the dm object 60 . nsubnet - local number of subnetworks 61 - Nsubnet - global number of subnetworks 62 63 Level: beginner 64 65 .seealso: DMNetworkCreate(), DMNetworkGetNumSubNetworks() 66 @*/ 67 PetscErrorCode DMNetworkSetNumSubNetworks(DM dm,PetscInt nsubnet,PetscInt Nsubnet) 68 { 69 PetscErrorCode ierr; 70 DM_Network *network = (DM_Network*)dm->data; 71 72 PetscFunctionBegin; 73 if (network->Nsubnet != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network"); 74 75 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 76 PetscValidLogicalCollectiveInt(dm,nsubnet,2); 77 PetscValidLogicalCollectiveInt(dm,Nsubnet,3); 78 79 if (Nsubnet == PETSC_DECIDE) { 80 if (nsubnet < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of local subnetworks %D cannot be less than 0",nsubnet); 81 ierr = MPIU_Allreduce(&nsubnet,&Nsubnet,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 82 } 83 if (Nsubnet < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Number of global subnetworks %D cannot be less than 1",Nsubnet); 84 85 network->Nsubnet = Nsubnet; 86 network->nsubnet = 0; /* initia value; will be determind by DMNetworkAddSubnetwork() */ 87 ierr = PetscCalloc1(Nsubnet,&network->subnet);CHKERRQ(ierr); 88 89 /* num of shared vertices */ 90 network->nsvtx = 0; 91 network->Nsvtx = 0; 92 PetscFunctionReturn(0); 93 } 94 95 /*@ 96 DMNetworkAddSubnetwork - Add a subnetwork 97 98 Collective on dm 99 100 Input Parameters: 101 + dm - the dm object 102 . name - name of the subnetwork 103 . nv - number of local vertices of this subnetwork 104 . ne - number of local edges of this subnetwork 105 - edgelist - list of edges for this subnetwork 106 107 Output Parameters: 108 . netnum - global index of the subnetwork 109 110 Notes: 111 There is no copy involved in this operation, only the pointer is referenced. The edgelist should 112 not be destroyed before the call to DMNetworkLayoutSetUp() 113 114 Level: beginner 115 116 Example usage: 117 Consider the following network: 118 .vb 119 network 1: v1 -> v2 -> v0 120 .ve 121 122 The resulting input 123 edgelist = [1 2 | 2 0] 124 125 .seealso: DMNetworkCreate(), DMNetworkSetNumSubnetworks() 126 @*/ 127 PetscErrorCode DMNetworkAddSubnetwork(DM dm,const char* name,PetscInt nv,PetscInt ne,PetscInt edgelist[],PetscInt *netnum) 128 { 129 PetscErrorCode ierr; 130 DM_Network *network = (DM_Network*)dm->data; 131 PetscInt i = network->nsubnet,a[2],b[2]; 132 133 PetscFunctionBegin; 134 if (name) { 135 ierr = PetscStrcpy(network->subnet[i].name,name);CHKERRQ(ierr); 136 } 137 138 network->subnet[i].nvtx = nv; 139 network->subnet[i].nedge = ne; 140 network->subnet[i].edgelist = edgelist; 141 142 /* Get global number of vertices and edges for subnet[i] */ 143 a[0] = nv; a[1] = ne; 144 ierr = MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 145 network->subnet[i].Nvtx = b[0]; 146 network->subnet[i].Nedge = b[1]; 147 148 /* ---------------------------------------------------------- 149 p=v or e; 150 subnet[0].pStart = 0 151 subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i]) 152 ----------------------------------------------------------------------- */ 153 /* GLOBAL subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */ 154 network->subnet[i].vStart = network->NVertices; 155 network->subnet[i].vEnd = network->subnet[i].vStart + network->subnet[i].Nvtx; /* global vEnd of subnet[i] */ 156 157 network->nVertices += nv; 158 network->NVertices += network->subnet[i].Nvtx; 159 160 /* LOCAL subnet[].eStart and eEnd, used by DMNetworkLayoutSetUp() */ 161 network->subnet[i].eStart = network->nEdges; 162 network->subnet[i].eEnd = network->subnet[i].eStart + ne; 163 network->nEdges += ne; 164 network->NEdges += network->subnet[i].Nedge; 165 166 ierr = PetscStrcpy(network->subnet[i].name,name);CHKERRQ(ierr); 167 if (netnum) *netnum = network->nsubnet; 168 network->nsubnet++; 169 PetscFunctionReturn(0); 170 } 171 172 /* 173 SetUp a single svtx struct. See SVtx defined in dmnetworkimpl.h 174 Set gidx and type if input v=(net,idx) is a from_vertex; 175 Get gid, type and index in the svtx array if input v=(net,idx) is a to_vertex. 176 177 Input: Nsvtx, svtx, net, idx, gidx 178 Output: gidx, svtype, svtx_idx 179 */ 180 static PetscErrorCode SVtxSetUp(PetscInt Nsvtx,SVtx *svtx,PetscInt net,PetscInt idx,PetscInt *gidx,SVtxType *svtype,PetscInt *svtx_idx) 181 { 182 PetscInt i,j,*svto; 183 SVtxType vtype; 184 185 PetscFunctionBegin; 186 if (!Nsvtx) PetscFunctionReturn(0); 187 188 vtype = SVNONE; 189 for (i=0; i<Nsvtx; i++) { 190 if (net == svtx[i].sv[0] && idx == svtx[i].sv[1]) { 191 /* (1) input vertex net.idx is a shared from_vertex, set its global index and output its svtype */ 192 svtx[i].gidx = *gidx; /* set gidx */ 193 vtype = SVFROM; 194 } else { /* loop over svtx[i].n */ 195 for (j=1; j<svtx[i].n; j++) { 196 svto = svtx[i].sv + 2*j; 197 if (net == svto[0] && idx == svto[1]) { 198 /* input vertex net.idx is a shared to_vertex, output its global index and its svtype */ 199 *gidx = svtx[i].gidx; /* output gidx for to_vertex */ 200 vtype = SVTO; 201 } 202 } 203 } 204 if (vtype != SVNONE) break; 205 } 206 if (svtype) *svtype = vtype; 207 if (svtx_idx) *svtx_idx = i; 208 PetscFunctionReturn(0); 209 } 210 211 /* 212 Add a new shared vertex sv=(net,idx) to table svtas[ita] 213 */ 214 static PetscErrorCode TableAddSVtx(PetscTable *svtas,PetscInt ita,PetscInt* tdata,PetscInt *sv_wk,PetscInt *ii,PetscInt *sedgelist,PetscInt k,DM_Network *network,PetscInt **ta2sv) 215 { 216 PetscInt net,idx,gidx,i=*ii; 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 net = sv_wk[2*i] = sedgelist[k]; 221 idx = sv_wk[2*i+1] = sedgelist[k+1]; 222 gidx = network->subnet[net].vStart + idx; 223 ierr = PetscTableAdd(svtas[ita],gidx+1,tdata[ita]+1,INSERT_VALUES);CHKERRQ(ierr); 224 *(ta2sv[ita] + tdata[ita]) = i; /* maps tdata to index of sv_wk; sv_wk keeps (net,idx) info */ 225 tdata[ita]++; (*ii)++; 226 PetscFunctionReturn(0); 227 } 228 229 /* 230 Create an array of shared vertices. See SVtx and SVtxType in dmnetworkimpl.h 231 232 Input: dm, Nsedgelist, sedgelist 233 Output: Nsvtx,svtx 234 235 Note: Output svtx is organized as 236 sv(net[0],idx[0]) --> sv(net[1],idx[1]) 237 --> sv(net[1],idx[1]) 238 ... 239 --> sv(net[n-1],idx[n-1]) 240 and net[0] < net[1] < ... < net[n-1] 241 where sv[0] has SVFROM type, sv[i], i>0, has SVTO type. 242 */ 243 static PetscErrorCode SVtxCreate(DM dm,PetscInt Nsedgelist,PetscInt *sedgelist,PetscInt *Nsvtx,SVtx **svtx) 244 { 245 PetscErrorCode ierr; 246 SVtx *sedges = NULL; 247 PetscInt *sv,k,j,nsv,*tdata,**ta2sv; 248 PetscTable *svtas; 249 PetscInt gidx,net,idx,i,nta,ita,idx_from,idx_to,n,*sv_wk; 250 DM_Network *network = (DM_Network*)dm->data; 251 PetscTablePosition ppos; 252 253 PetscFunctionBegin; 254 /* (1) Crete ctables svtas */ 255 ierr = PetscCalloc4(Nsedgelist,&svtas,Nsedgelist,&tdata,4*Nsedgelist,&sv_wk,2*Nsedgelist,&ta2sv);CHKERRQ(ierr); 256 257 j = 0; /* sedgelist counter */ 258 k = 0; /* sedgelist vertex counter j = 4*k */ 259 i = 0; /* sv_wk (vertices added to the ctables) counter */ 260 nta = 0; /* num of sv tables created */ 261 262 /* for j=0 */ 263 ierr = PetscTableCreate(2*Nsedgelist,network->NVertices+1,&svtas[nta]);CHKERRQ(ierr); 264 ierr = PetscMalloc1(2*Nsedgelist,&ta2sv[nta]);CHKERRQ(ierr); 265 266 ierr = TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);CHKERRQ(ierr); 267 ierr = TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);CHKERRQ(ierr); 268 nta++; k += 4; 269 270 for (j = 1; j < Nsedgelist; j++) { 271 for (ita = 0; ita < nta; ita++) { 272 /* vfrom */ 273 net = sedgelist[k]; idx = sedgelist[k+1]; 274 gidx = network->subnet[net].vStart + idx; /* global index of the vertex net.idx before merging shared vertices */ 275 ierr = PetscTableFind(svtas[ita],gidx+1,&idx_from);CHKERRQ(ierr); 276 277 /* vto */ 278 net = sedgelist[k+2]; idx = sedgelist[k+3]; 279 gidx = network->subnet[net].vStart + idx; 280 ierr = PetscTableFind(svtas[ita],gidx+1,&idx_to);CHKERRQ(ierr); 281 282 if (idx_from || idx_to) { /* vfrom or vto is on table svtas[ita] */ 283 idx_from--; idx_to--; 284 if (idx_from < 0) { /* vto is on svtas[ita] */ 285 ierr = TableAddSVtx(svtas,ita,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);CHKERRQ(ierr); 286 break; 287 } else if (idx_to < 0) { 288 ierr = TableAddSVtx(svtas,ita,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);CHKERRQ(ierr); 289 break; 290 } 291 } 292 } 293 294 if (ita == nta) { 295 ierr = PetscTableCreate(2*Nsedgelist,network->NVertices+1,&svtas[nta]);CHKERRQ(ierr); 296 ierr = PetscMalloc1(2*Nsedgelist, &ta2sv[nta]);CHKERRQ(ierr); 297 298 ierr = TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);CHKERRQ(ierr); 299 ierr = TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);CHKERRQ(ierr); 300 nta++; 301 } 302 k += 4; 303 } 304 305 /* (2) Construct sedges from ctable 306 sedges: edges connect vertex sv[0]=(net[0],idx[0]) to vertices sv[k], k=1,...,n-1; 307 net[k], k=0, ...,n-1, are in ascending order */ 308 ierr = PetscMalloc1(nta,&sedges);CHKERRQ(ierr); 309 for (nsv = 0; nsv < nta; nsv++) { 310 /* for a single svtx, put shared vertices in ascending order of gidx */ 311 ierr = PetscTableGetCount(svtas[nsv],&n);CHKERRQ(ierr); 312 ierr = PetscCalloc1(2*n,&sv);CHKERRQ(ierr); 313 sedges[nsv].sv = sv; 314 sedges[nsv].n = n; 315 sedges[nsv].gidx = -1; /* initialization */ 316 317 ierr = PetscTableGetHeadPosition(svtas[nsv],&ppos);CHKERRQ(ierr); 318 for (k=0; k<n; k++) { /* gidx is sorted in ascending order */ 319 ierr = PetscTableGetNext(svtas[nsv],&ppos,&gidx,&i);CHKERRQ(ierr); 320 gidx--; i--; 321 322 j = ta2sv[nsv][i]; /* maps i to index of sv_wk */ 323 sv[2*k] = sv_wk[2*j]; 324 sv[2*k+1] = sv_wk[2*j + 1]; 325 } 326 } 327 328 for (j=0; j<nta; j++) { 329 ierr = PetscTableDestroy(&svtas[j]);CHKERRQ(ierr); 330 ierr = PetscFree(ta2sv[j]);CHKERRQ(ierr); 331 } 332 ierr = PetscFree4(svtas,tdata,sv_wk,ta2sv);CHKERRQ(ierr); 333 334 *Nsvtx = nta; 335 *svtx = sedges; 336 PetscFunctionReturn(0); 337 } 338 339 static PetscErrorCode DMNetworkLayoutSetUp_Coupling(DM dm) 340 { 341 PetscErrorCode ierr; 342 DM_Network *network = (DM_Network*)dm->data; 343 PetscInt i,j,ctr,np,*edges,*subnetvtx,vStart; 344 PetscInt k,*vidxlTog,Nsv=0,Nsubnet=network->Nsubnet; 345 PetscInt *sedgelist=network->sedgelist; 346 const PetscInt *cone; 347 MPI_Comm comm; 348 PetscMPIInt size,rank,*recvcounts=NULL,*displs=NULL; 349 PetscInt net,idx,gidx,nmerged,e,v,vfrom,vto,*vrange,*eowners,gidx_from,net_from,sv_idx; 350 SVtxType svtype = SVNONE; 351 SVtx *svtx=NULL; 352 PetscSection sectiong; 353 354 PetscFunctionBegin; 355 /* This implementation requires user input each subnet by a single processor, thus subnet[net].nvtx=subnet[net].Nvtx */ 356 for (net=0; net<Nsubnet; net++) { 357 if (network->subnet[net].nvtx && network->subnet[net].nvtx != network->subnet[net].Nvtx) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"subnetwork %D local num of vertices %D != %D global num",net,network->subnet[net].nvtx,network->subnet[net].Nvtx); 358 } 359 360 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 361 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 362 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 363 364 /* (1) Create svtx[] from sedgelist */ 365 /* -------------------------------- */ 366 /* Nsv: global number of SVtx; svtx: shared vertices, see SVtx in dmnetworkimpl.h */ 367 ierr = SVtxCreate(dm,network->Nsvtx,sedgelist,&Nsv,&svtx);CHKERRQ(ierr); 368 369 /* (2) Setup svtx; Shared vto vertices are merged to their vfrom vertex with same global vetex index (gidx) */ 370 /* -------------------------------------------------------------------------------------------------------- */ 371 /* (2.1) compute vrage[rank]: global index of 1st local vertex in proc[rank] */ 372 ierr = PetscMalloc3(size+1,&vrange,size,&displs,size,&recvcounts);CHKERRQ(ierr); 373 for (i=0; i<size; i++) {displs[i] = i; recvcounts[i] = 1;} 374 375 vrange[0] = 0; 376 ierr = MPI_Allgatherv(&network->nVertices,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);CHKERRMPI(ierr); 377 for (i=2; i<size+1; i++) { 378 vrange[i] += vrange[i-1]; 379 } 380 381 /* (2.2) Create vidxlTog: maps UN-MERGED local vertex index i to global index gidx (plex, excluding ghost vertices) */ 382 ierr = PetscMalloc1(network->nVertices,&vidxlTog);CHKERRQ(ierr); 383 i = 0; gidx = 0; 384 nmerged = 0; /* local num of merged vertices */ 385 network->nsvtx = 0; 386 for (net=0; net<Nsubnet; net++) { 387 for (idx=0; idx<network->subnet[net].Nvtx; idx++) { 388 gidx_from = gidx; 389 sv_idx = -1; 390 391 ierr = SVtxSetUp(Nsv,svtx,net,idx,&gidx_from,&svtype,&sv_idx);CHKERRQ(ierr); 392 if (svtype == SVTO) { 393 if (network->subnet[net].nvtx) {/* this proc owns sv_to */ 394 net_from = svtx[sv_idx].sv[0]; /* subnet num of its shared vertex */ 395 if (network->subnet[net_from].nvtx == 0) { 396 /* this proc does not own v_from, thus a new local coupling vertex */ 397 network->nsvtx++; 398 } 399 vidxlTog[i++] = gidx_from; 400 nmerged++; /* a coupling vertex -- merged */ 401 } 402 } else { 403 if (svtype == SVFROM) { 404 if (network->subnet[net].nvtx) { 405 /* this proc owns this v_from, a new local coupling vertex */ 406 network->nsvtx++; 407 } 408 } 409 if (network->subnet[net].nvtx) vidxlTog[i++] = gidx; 410 gidx++; 411 } 412 } 413 } 414 #if defined(PETSC_USE_DEBUG) 415 if (i != network->nVertices) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"%D != %D nVertices",i,network->nVertices); 416 #endif 417 418 /* (2.3) Setup svtable for querry shared vertices */ 419 for (v=0; v<Nsv; v++) { 420 gidx = svtx[v].gidx; 421 ierr = PetscTableAdd(network->svtable,gidx+1,v+1,INSERT_VALUES);CHKERRQ(ierr); 422 } 423 424 /* (2.4) Shared vertices in the subnetworks are merged, update global NVertices: np = sum(local nmerged) */ 425 ierr = MPI_Allreduce(&nmerged,&np,1,MPIU_INT,MPI_SUM,comm);CHKERRMPI(ierr); 426 network->NVertices -= np; 427 428 ierr = PetscCalloc1(2*network->nEdges,&edges);CHKERRQ(ierr); 429 ierr = PetscCalloc1(network->nVertices+network->nsvtx,&network->subnetvtx);CHKERRQ(ierr); 430 431 ctr = 0; 432 for (net=0; net<Nsubnet; net++) { 433 for (j = 0; j < network->subnet[net].nedge; j++) { 434 /* vfrom: */ 435 i = network->subnet[net].edgelist[2*j] + (network->subnet[net].vStart - vrange[rank]); 436 edges[2*ctr] = vidxlTog[i]; 437 438 /* vto */ 439 i = network->subnet[net].edgelist[2*j+1] + (network->subnet[net].vStart - vrange[rank]); 440 edges[2*ctr+1] = vidxlTog[i]; 441 ctr++; 442 } 443 } 444 ierr = PetscFree3(vrange,displs,recvcounts);CHKERRQ(ierr); 445 ierr = PetscFree(vidxlTog);CHKERRQ(ierr); 446 447 /* (3) Create network->plex */ 448 ierr = DMCreate(comm,&network->plex);CHKERRQ(ierr); 449 ierr = DMSetType(network->plex,DMPLEX);CHKERRQ(ierr); 450 ierr = DMSetDimension(network->plex,1);CHKERRQ(ierr); 451 if (size == 1) { 452 ierr = DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices-nmerged,2,edges);CHKERRQ(ierr); 453 } else { 454 ierr = DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices-nmerged,PETSC_DECIDE,2,edges,NULL);CHKERRQ(ierr); 455 } 456 457 ierr = DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);CHKERRQ(ierr); 458 ierr = DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);CHKERRQ(ierr); 459 ierr = DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);CHKERRQ(ierr); 460 vStart = network->vStart; 461 462 ierr = PetscSectionCreate(comm,&network->DataSection);CHKERRQ(ierr); 463 ierr = PetscSectionCreate(comm,&network->DofSection);CHKERRQ(ierr); 464 ierr = PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);CHKERRQ(ierr); 465 ierr = PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);CHKERRQ(ierr); 466 467 network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType); 468 np = network->pEnd - network->pStart; 469 ierr = PetscCalloc2(np,&network->header,np,&network->cvalue);CHKERRQ(ierr); 470 471 /* (4) Create vidxlTog: maps MERGED plex local vertex index (including ghosts) to User's global vertex index (without merging shared vertices) */ 472 np = network->vEnd - vStart; /* include ghost vertices */ 473 ierr = PetscMalloc2(np,&vidxlTog,size+1,&eowners);CHKERRQ(ierr); 474 475 ctr = 0; 476 for (e=network->eStart; e<network->eEnd; e++) { 477 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 478 vidxlTog[cone[0] - vStart] = edges[2*ctr]; 479 vidxlTog[cone[1] - vStart] = edges[2*ctr+1]; 480 ctr++; 481 } 482 ierr = PetscFree(edges);CHKERRQ(ierr); 483 484 /* (5) Create vertices and edges array for the subnetworks */ 485 subnetvtx = network->subnetvtx; 486 for (j=0; j < Nsubnet; j++) { 487 ierr = PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);CHKERRQ(ierr); 488 network->subnet[j].vertices = subnetvtx; 489 subnetvtx += network->subnet[j].nvtx; 490 } 491 network->svertices = subnetvtx; 492 493 /* Get edge ownership */ 494 np = network->eEnd - network->eStart; /* num of local edges */ 495 ierr = MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);CHKERRMPI(ierr); 496 eowners[0] = 0; 497 for (i=2; i<=size; i++) eowners[i] += eowners[i-1]; 498 499 e = 0; 500 for (i=0; i < Nsubnet; i++) { 501 v = 0; 502 for (j = 0; j < network->subnet[i].nedge; j++) { 503 504 /* edge e */ 505 network->header[e].index = e + eowners[rank]; /* Global edge index */ 506 network->header[e].subnetid = i; /* Subnetwork id */ 507 network->subnet[i].edges[j] = e; 508 network->header[e].ndata = 0; 509 network->header[e].offset[0] = 0; 510 network->header[e].offsetvarrel[0] = 0; 511 ierr = PetscSectionAddDof(network->DataSection,e,network->dataheadersize);CHKERRQ(ierr); 512 513 /* connected vertices */ 514 ierr = DMPlexGetCone(network->plex,e,&cone);CHKERRQ(ierr); 515 516 /* vertex cone[0] */ 517 vfrom = network->subnet[i].edgelist[2*v]; 518 network->header[cone[0]].index = vidxlTog[cone[0]-vStart]; /* Global vertex index */ 519 network->header[cone[0]].subnetid = i; /* Subnetwork id */ 520 network->subnet[i].vertices[vfrom] = cone[0]; /* user's subnet[].dix = petsc's v */ 521 522 /* vertex cone[1] */ 523 vto = network->subnet[i].edgelist[2*v+1]; 524 network->header[cone[1]].index = vidxlTog[cone[1]-vStart]; /* Global vertex index */ 525 network->header[cone[1]].subnetid = i; 526 network->subnet[i].vertices[vto]= cone[1]; 527 528 e++; v++; 529 } 530 } 531 532 /* Set vertex array for the subnetworks */ 533 k = 0; 534 for (v=vStart; v<network->vEnd; v++) { /* local vertices, including ghosts */ 535 network->header[v].ndata = 0; 536 network->header[v].offset[0] = 0; 537 network->header[v].offsetvarrel[0] = 0; 538 ierr = PetscSectionAddDof(network->DataSection,v,network->dataheadersize);CHKERRQ(ierr); 539 540 /* shared vertex */ 541 ierr = PetscTableFind(network->svtable,vidxlTog[v-vStart]+1,&i);CHKERRQ(ierr); 542 if (i) network->svertices[k++] = v; 543 } 544 #if defined(PETSC_USE_DEBUG) 545 if (k != network->nsvtx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"k %D != %D nsvtx",k,network->nsvtx); 546 #endif 547 548 ierr = PetscFree2(vidxlTog,eowners);CHKERRQ(ierr); 549 550 network->svtx = svtx; 551 network->Nsvtx = Nsv; 552 ierr = PetscFree(sedgelist);CHKERRQ(ierr); 553 554 /* Create a global section to be used by DMNetworkIsGhostVertex() which is a non-collective routine */ 555 ierr = DMGetGlobalSection(network->plex,§iong);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 /*@ 560 DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network 561 562 Collective on dm 563 564 Input Parameters: 565 . DM - the dmnetwork object 566 567 Notes: 568 This routine should be called after the network sizes and edgelists have been provided. It creates 569 the bare layout of the network and sets up the network to begin insertion of components. 570 571 All the components should be registered before calling this routine. 572 573 Level: beginner 574 575 .seealso: DMNetworkSetNumSubNetworks(), DMNetworkAddSubnetwork() 576 @*/ 577 PetscErrorCode DMNetworkLayoutSetUp(DM dm) 578 { 579 PetscErrorCode ierr; 580 DM_Network *network = (DM_Network*)dm->data; 581 PetscInt i,j,ctr,Nsubnet=network->Nsubnet,*eowners,np,*edges,*subnetvtx; 582 PetscInt e,v,vfrom,vto; 583 const PetscInt *cone; 584 MPI_Comm comm; 585 PetscMPIInt size,rank; 586 587 PetscFunctionBegin; 588 if (network->nsubnet != network->Nsubnet) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Must call DMNetworkAddSubnetwork() %D times",network->Nsubnet); 589 590 /* Create svtable for querry shared vertices */ 591 ierr = PetscTableCreate(network->Nsvtx,network->NVertices+1,&network->svtable);CHKERRQ(ierr); 592 593 if (network->Nsvtx) { 594 ierr = DMNetworkLayoutSetUp_Coupling(dm);CHKERRQ(ierr); 595 PetscFunctionReturn(0); 596 } 597 598 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 599 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 600 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 601 602 /* Create LOCAL edgelist for the network by concatenating local input edgelists of the subnetworks */ 603 ierr = PetscCalloc1(2*network->nEdges,&edges);CHKERRQ(ierr); 604 ctr = 0; 605 for (i=0; i < Nsubnet; i++) { 606 for (j = 0; j < network->subnet[i].nedge; j++) { 607 edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j]; 608 edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1]; 609 ctr++; 610 } 611 } 612 613 /* Create network->plex; One dimensional network, numCorners=2 */ 614 ierr = DMCreate(comm,&network->plex);CHKERRQ(ierr); 615 ierr = DMSetType(network->plex,DMPLEX);CHKERRQ(ierr); 616 ierr = DMSetDimension(network->plex,1);CHKERRQ(ierr); 617 if (size == 1) { 618 ierr = DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices,2,edges);CHKERRQ(ierr); 619 } else { 620 ierr = DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices,PETSC_DECIDE,2,edges,NULL);CHKERRQ(ierr); 621 } 622 ierr = PetscFree(edges);CHKERRQ(ierr); /* local edge list with global idx used by DMPlexBuildFromCellList() */ 623 624 ierr = DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);CHKERRQ(ierr); 625 ierr = DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);CHKERRQ(ierr); 626 ierr = DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);CHKERRQ(ierr); 627 628 ierr = PetscSectionCreate(comm,&network->DataSection);CHKERRQ(ierr); 629 ierr = PetscSectionCreate(comm,&network->DofSection);CHKERRQ(ierr); 630 ierr = PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);CHKERRQ(ierr); 631 ierr = PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);CHKERRQ(ierr); 632 633 network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType); 634 np = network->pEnd - network->pStart; 635 ierr = PetscCalloc2(np,&network->header,np,&network->cvalue);CHKERRQ(ierr); 636 637 /* Create edge and vertex arrays for the subnetworks */ 638 for (j=0; j < network->Nsubnet; j++) { 639 ierr = PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);CHKERRQ(ierr); 640 } 641 642 /* Get edge ownership */ 643 ierr = PetscMalloc1(size+1,&eowners);CHKERRQ(ierr); 644 np = network->eEnd - network->eStart; 645 ierr = MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);CHKERRMPI(ierr); 646 eowners[0] = 0; 647 for (i=2; i<=size; i++) eowners[i] += eowners[i-1]; 648 649 /* Set network->subnet[*].vertices on array network->subnetvtx */ 650 np = 0; 651 for (j=0; j<network->Nsubnet; j++) { 652 /* sum up subnet[j].Nvtx instead of subnet[j].nvtx, because a subnet might be owned by more than one processor; 653 below, subnet[i].vertices[vfrom/vto] requires vfrom/vto =0, ...,Nvtx-1 654 */ 655 if (network->subnet[j].nvtx) np += network->subnet[j].Nvtx; 656 } 657 658 ierr = PetscCalloc1(np,&network->subnetvtx);CHKERRQ(ierr); /* Maps local vertex to local subnetwork's vertex */ 659 subnetvtx = network->subnetvtx; 660 for (j=0; j<network->Nsubnet; j++) { 661 network->subnet[j].vertices = subnetvtx; 662 if (network->subnet[j].nvtx) subnetvtx += network->subnet[j].Nvtx; 663 } 664 665 /* Setup edge and vertex arrays for subnetworks */ 666 e = 0; 667 for (i=0; i < Nsubnet; i++) { 668 v = 0; 669 for (j = 0; j < network->subnet[i].nedge; j++) { 670 /* edge e */ 671 network->header[e].index = e + eowners[rank]; /* Global edge index */ 672 network->header[e].subnetid = i; 673 network->subnet[i].edges[j] = e; 674 675 network->header[e].ndata = 0; 676 network->header[e].offset[0] = 0; 677 network->header[e].offsetvarrel[0] = 0; 678 ierr = PetscSectionAddDof(network->DataSection,e,network->dataheadersize);CHKERRQ(ierr); 679 680 /* connected vertices */ 681 ierr = DMPlexGetCone(network->plex,e,&cone);CHKERRQ(ierr); 682 683 /* vertex cone[0] */ 684 vfrom = network->subnet[i].edgelist[2*v]; /* =subnet[i].idx, Global index! */ 685 network->header[cone[0]].index = vfrom + network->subnet[i].vStart; /* Global vertex index */ 686 network->header[cone[0]].subnetid = i; /* Subnetwork id */ 687 network->subnet[i].vertices[vfrom] = cone[0]; /* user's subnet[].dix = petsc's v */ 688 689 /* vertex cone[1] */ 690 vto = network->subnet[i].edgelist[2*v+1]; /* =subnet[i].idx, Global index! */ 691 network->header[cone[1]].index = vto + network->subnet[i].vStart; /* Global vertex index */ 692 network->header[cone[1]].subnetid = i; 693 network->subnet[i].vertices[vto] = cone[1]; /* user's subnet[].dix = petsc's v */ 694 695 e++; v++; 696 } 697 } 698 ierr = PetscFree(eowners);CHKERRQ(ierr); 699 700 for (v = network->vStart; v < network->vEnd; v++) { 701 network->header[v].ndata = 0; 702 network->header[v].offset[0] = 0; 703 network->header[v].offsetvarrel[0] = 0; 704 ierr = PetscSectionAddDof(network->DataSection,v,network->dataheadersize);CHKERRQ(ierr); 705 } 706 PetscFunctionReturn(0); 707 } 708 709 /*@C 710 DMNetworkGetSubnetwork - Returns the information about a requested subnetwork 711 712 Not collective 713 714 Input Parameters: 715 + dm - the DM object 716 - netnum - the global index of the subnetwork 717 718 Output Parameters: 719 + nv - number of vertices (local) 720 . ne - number of edges (local) 721 . vtx - local vertices of the subnetwork 722 - edge - local edges of the subnetwork 723 724 Notes: 725 Cannot call this routine before DMNetworkLayoutSetup() 726 727 Level: intermediate 728 729 .seealso: DMNetworkCreate(), DMNetworkAddSubnetwork(), DMNetworkLayoutSetUp() 730 @*/ 731 PetscErrorCode DMNetworkGetSubnetwork(DM dm,PetscInt netnum,PetscInt *nv,PetscInt *ne,const PetscInt **vtx,const PetscInt **edge) 732 { 733 DM_Network *network = (DM_Network*)dm->data; 734 735 PetscFunctionBegin; 736 if (netnum >= network->Nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet index %D exceeds the num of subnets %D",netnum,network->Nsubnet); 737 if (nv) *nv = network->subnet[netnum].nvtx; 738 if (ne) *ne = network->subnet[netnum].nedge; 739 if (vtx) *vtx = network->subnet[netnum].vertices; 740 if (edge) *edge = network->subnet[netnum].edges; 741 PetscFunctionReturn(0); 742 } 743 744 /*@ 745 DMNetworkAddSharedVertices - Add shared vertices that connect two given subnetworks 746 747 Collective on dm 748 749 Input Parameters: 750 + dm - the dm object 751 . anetnum - first subnetwork global numbering returned by DMNetworkAddSubnetwork() 752 . bnetnum - second subnetwork global numbering returned by DMNetworkAddSubnetwork() 753 . nsvtx - number of vertices that are shared by the two subnetworks 754 . asvtx - vertex index in the first subnetwork 755 - bsvtx - vertex index in the second subnetwork 756 757 Level: beginner 758 759 .seealso: DMNetworkCreate(), DMNetworkAddSubnetwork(), DMNetworkGetSharedVertices() 760 @*/ 761 PetscErrorCode DMNetworkAddSharedVertices(DM dm,PetscInt anetnum,PetscInt bnetnum,PetscInt nsvtx,PetscInt asvtx[],PetscInt bsvtx[]) 762 { 763 PetscErrorCode ierr; 764 DM_Network *network = (DM_Network*)dm->data; 765 PetscInt i,nsubnet = network->Nsubnet,*sedgelist,Nsvtx=network->Nsvtx; 766 767 PetscFunctionBegin; 768 if (anetnum == bnetnum) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Subnetworks must have different netnum"); 769 if (anetnum < 0 || bnetnum < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"netnum cannot be negative"); 770 if (!Nsvtx) { 771 /* allocate network->sedgelist to hold at most 2*nsubnet pairs of shared vertices */ 772 ierr = PetscMalloc1(2*4*nsubnet,&network->sedgelist);CHKERRQ(ierr); 773 } 774 775 sedgelist = network->sedgelist; 776 for (i=0; i<nsvtx; i++) { 777 sedgelist[4*Nsvtx] = anetnum; sedgelist[4*Nsvtx+1] = asvtx[i]; 778 sedgelist[4*Nsvtx+2] = bnetnum; sedgelist[4*Nsvtx+3] = bsvtx[i]; 779 Nsvtx++; 780 } 781 if (Nsvtx > 2*nsubnet) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"allocate more space for coupling edgelist"); 782 network->Nsvtx = Nsvtx; 783 PetscFunctionReturn(0); 784 } 785 786 /*@C 787 DMNetworkGetSharedVertices - Returns the info for the shared vertices 788 789 Not collective 790 791 Input Parameters: 792 . dm - the DM object 793 794 Output Parameters: 795 + nsv - number of local shared vertices 796 - svtx - local shared vertices 797 798 Notes: 799 Cannot call this routine before DMNetworkLayoutSetup() 800 801 Level: intermediate 802 803 .seealso: DMNetworkGetSubnetwork(), DMNetworkLayoutSetUp(), DMNetworkAddSharedVertices() 804 @*/ 805 PetscErrorCode DMNetworkGetSharedVertices(DM dm,PetscInt *nsv,const PetscInt **svtx) 806 { 807 DM_Network *net = (DM_Network*)dm->data; 808 809 PetscFunctionBegin; 810 if (net->Nsvtx) { 811 *nsv = net->nsvtx; 812 *svtx = net->svertices; 813 } else { 814 *nsv = 0; 815 *svtx = NULL; 816 } 817 PetscFunctionReturn(0); 818 } 819 820 /*@C 821 DMNetworkRegisterComponent - Registers the network component 822 823 Logically collective on dm 824 825 Input Parameters: 826 + dm - the network object 827 . name - the component name 828 - size - the storage size in bytes for this component data 829 830 Output Parameters: 831 . key - an integer key that defines the component 832 833 Notes 834 This routine should be called by all processors before calling DMNetworkLayoutSetup(). 835 836 Level: beginner 837 838 .seealso: DMNetworkCreate(), DMNetworkLayoutSetUp() 839 @*/ 840 PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key) 841 { 842 PetscErrorCode ierr; 843 DM_Network *network = (DM_Network*) dm->data; 844 DMNetworkComponent *component=&network->component[network->ncomponent]; 845 PetscBool flg=PETSC_FALSE; 846 PetscInt i; 847 848 PetscFunctionBegin; 849 for (i=0; i < network->ncomponent; i++) { 850 ierr = PetscStrcmp(component->name,name,&flg);CHKERRQ(ierr); 851 if (flg) { 852 *key = i; 853 PetscFunctionReturn(0); 854 } 855 } 856 if (network->ncomponent == MAX_COMPONENTS) { 857 SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS); 858 } 859 860 ierr = PetscStrcpy(component->name,name);CHKERRQ(ierr); 861 component->size = size/sizeof(DMNetworkComponentGenericDataType); 862 *key = network->ncomponent; 863 network->ncomponent++; 864 PetscFunctionReturn(0); 865 } 866 867 /*@ 868 DMNetworkGetVertexRange - Get the bounds [start, end) for the local vertices 869 870 Not Collective 871 872 Input Parameters: 873 . dm - the DMNetwork object 874 875 Output Parameters: 876 + vStart - the first vertex point 877 - vEnd - one beyond the last vertex point 878 879 Level: beginner 880 881 .seealso: DMNetworkGetEdgeRange() 882 @*/ 883 PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd) 884 { 885 DM_Network *network = (DM_Network*)dm->data; 886 887 PetscFunctionBegin; 888 if (vStart) *vStart = network->vStart; 889 if (vEnd) *vEnd = network->vEnd; 890 PetscFunctionReturn(0); 891 } 892 893 /*@ 894 DMNetworkGetEdgeRange - Get the bounds [start, end) for the local edges 895 896 Not Collective 897 898 Input Parameters: 899 . dm - the DMNetwork object 900 901 Output Parameters: 902 + eStart - The first edge point 903 - eEnd - One beyond the last edge point 904 905 Level: beginner 906 907 .seealso: DMNetworkGetVertexRange() 908 @*/ 909 PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd) 910 { 911 DM_Network *network = (DM_Network*)dm->data; 912 913 PetscFunctionBegin; 914 if (eStart) *eStart = network->eStart; 915 if (eEnd) *eEnd = network->eEnd; 916 PetscFunctionReturn(0); 917 } 918 919 static PetscErrorCode DMNetworkGetIndex(DM dm,PetscInt p,PetscInt *index) 920 { 921 PetscErrorCode ierr; 922 DM_Network *network = (DM_Network*)dm->data; 923 PetscInt offsetp; 924 DMNetworkComponentHeader header; 925 926 PetscFunctionBegin; 927 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 928 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 929 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp); 930 *index = header->index; 931 PetscFunctionReturn(0); 932 } 933 934 /*@ 935 DMNetworkGetGlobalEdgeIndex - Get the global numbering for the edge on the network 936 937 Not Collective 938 939 Input Parameters: 940 + dm - DMNetwork object 941 - p - edge point 942 943 Output Parameters: 944 . index - the global numbering for the edge 945 946 Level: intermediate 947 948 .seealso: DMNetworkGetGlobalVertexIndex() 949 @*/ 950 PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index) 951 { 952 PetscErrorCode ierr; 953 954 PetscFunctionBegin; 955 ierr = DMNetworkGetIndex(dm,p,index);CHKERRQ(ierr); 956 PetscFunctionReturn(0); 957 } 958 959 /*@ 960 DMNetworkGetGlobalVertexIndex - Get the global numbering for the vertex on the network 961 962 Not Collective 963 964 Input Parameters: 965 + dm - DMNetwork object 966 - p - vertex point 967 968 Output Parameters: 969 . index - the global numbering for the vertex 970 971 Level: intermediate 972 973 .seealso: DMNetworkGetGlobalEdgeIndex() 974 @*/ 975 PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index) 976 { 977 PetscErrorCode ierr; 978 979 PetscFunctionBegin; 980 ierr = DMNetworkGetIndex(dm,p,index);CHKERRQ(ierr); 981 PetscFunctionReturn(0); 982 } 983 984 /*@ 985 DMNetworkGetNumComponents - Get the number of components at a vertex/edge 986 987 Not Collective 988 989 Input Parameters: 990 + dm - the DMNetwork object 991 - p - vertex/edge point 992 993 Output Parameters: 994 . numcomponents - Number of components at the vertex/edge 995 996 Level: beginner 997 998 .seealso: DMNetworkRegisterComponent(), DMNetworkAddComponent() 999 @*/ 1000 PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents) 1001 { 1002 PetscErrorCode ierr; 1003 PetscInt offset; 1004 DM_Network *network = (DM_Network*)dm->data; 1005 1006 PetscFunctionBegin; 1007 ierr = PetscSectionGetOffset(network->DataSection,p,&offset);CHKERRQ(ierr); 1008 *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata; 1009 PetscFunctionReturn(0); 1010 } 1011 1012 /*@ 1013 DMNetworkGetLocalVecOffset - Get the offset for accessing the variables associated with a component at the given vertex/edge from the local vector 1014 1015 Not Collective 1016 1017 Input Parameters: 1018 + dm - the DMNetwork object 1019 . p - the edge/vertex point 1020 - compnum - component number; use ALL_COMPONENTS if no specific component is requested 1021 1022 Output Parameters: 1023 . offset - the local offset 1024 1025 Level: intermediate 1026 1027 .seealso: DMGetLocalVector(), DMNetworkGetComponent(), DMNetworkGetGlobalVecOffset() 1028 @*/ 1029 PetscErrorCode DMNetworkGetLocalVecOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset) 1030 { 1031 PetscErrorCode ierr; 1032 DM_Network *network = (DM_Network*)dm->data; 1033 PetscInt offsetp,offsetd; 1034 DMNetworkComponentHeader header; 1035 1036 PetscFunctionBegin; 1037 ierr = PetscSectionGetOffset(network->plex->localSection,p,&offsetp);CHKERRQ(ierr); 1038 if (compnum == ALL_COMPONENTS) { 1039 *offset = offsetp; 1040 PetscFunctionReturn(0); 1041 } 1042 1043 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetd);CHKERRQ(ierr); 1044 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd); 1045 *offset = offsetp + header->offsetvarrel[compnum]; 1046 PetscFunctionReturn(0); 1047 } 1048 1049 /*@ 1050 DMNetworkGetGlobalVecOffset - Get the global offset for accessing the variables associated with a component for the given vertex/edge from the global vector 1051 1052 Not Collective 1053 1054 Input Parameters: 1055 + dm - the DMNetwork object 1056 . p - the edge/vertex point 1057 - compnum - component number; use ALL_COMPONENTS if no specific component is requested 1058 1059 Output Parameters: 1060 . offsetg - the global offset 1061 1062 Level: intermediate 1063 1064 .seealso: DMNetworkGetLocalVecOffset(), DMGetGlobalVector(), DMNetworkGetComponent() 1065 @*/ 1066 PetscErrorCode DMNetworkGetGlobalVecOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg) 1067 { 1068 PetscErrorCode ierr; 1069 DM_Network *network = (DM_Network*)dm->data; 1070 PetscInt offsetp,offsetd; 1071 DMNetworkComponentHeader header; 1072 1073 PetscFunctionBegin; 1074 ierr = PetscSectionGetOffset(network->plex->globalSection,p,&offsetp);CHKERRQ(ierr); 1075 if (offsetp < 0) offsetp = -(offsetp + 1); /* Convert to actual global offset for ghost vertex */ 1076 1077 if (compnum == ALL_COMPONENTS) { 1078 *offsetg = offsetp; 1079 PetscFunctionReturn(0); 1080 } 1081 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetd);CHKERRQ(ierr); 1082 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd); 1083 *offsetg = offsetp + header->offsetvarrel[compnum]; 1084 PetscFunctionReturn(0); 1085 } 1086 1087 /*@ 1088 DMNetworkGetEdgeOffset - Get the offset for accessing the variables associated with the given edge from the local subvector 1089 1090 Not Collective 1091 1092 Input Parameters: 1093 + dm - the DMNetwork object 1094 - p - the edge point 1095 1096 Output Parameters: 1097 . offset - the offset 1098 1099 Level: intermediate 1100 1101 .seealso: DMNetworkGetLocalVecOffset(), DMGetLocalVector() 1102 @*/ 1103 PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset) 1104 { 1105 PetscErrorCode ierr; 1106 DM_Network *network = (DM_Network*)dm->data; 1107 1108 PetscFunctionBegin; 1109 ierr = PetscSectionGetOffset(network->edge.DofSection,p,offset);CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 1113 /*@ 1114 DMNetworkGetVertexOffset - Get the offset for accessing the variables associated with the given vertex from the local subvector 1115 1116 Not Collective 1117 1118 Input Parameters: 1119 + dm - the DMNetwork object 1120 - p - the vertex point 1121 1122 Output Parameters: 1123 . offset - the offset 1124 1125 Level: intermediate 1126 1127 .seealso: DMNetworkGetEdgeOffset(), DMGetLocalVector() 1128 @*/ 1129 PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset) 1130 { 1131 PetscErrorCode ierr; 1132 DM_Network *network = (DM_Network*)dm->data; 1133 1134 PetscFunctionBegin; 1135 p -= network->vStart; 1136 ierr = PetscSectionGetOffset(network->vertex.DofSection,p,offset);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /*@ 1141 DMNetworkAddComponent - Adds a network component and number of variables at the given point (vertex/edge) 1142 1143 Not Collective 1144 1145 Input Parameters: 1146 + dm - the DMNetworkObject 1147 . p - the vertex/edge point 1148 . componentkey - component key returned while registering the component; ignored if compvalue=NULL 1149 . compvalue - pointer to the data structure for the component, or NULL if not required. 1150 - nvar - number of variables for the component at the vertex/edge point 1151 1152 Level: beginner 1153 1154 .seealso: DMNetworkGetComponent() 1155 @*/ 1156 PetscErrorCode DMNetworkAddComponent(DM dm,PetscInt p,PetscInt componentkey,void* compvalue,PetscInt nvar) 1157 { 1158 PetscErrorCode ierr; 1159 DM_Network *network = (DM_Network*)dm->data; 1160 DMNetworkComponent *component = &network->component[componentkey]; 1161 DMNetworkComponentHeader header = &network->header[p]; 1162 DMNetworkComponentValue cvalue = &network->cvalue[p]; 1163 PetscBool sharedv=PETSC_FALSE; 1164 PetscInt compnum=header->ndata; 1165 1166 PetscFunctionBegin; 1167 ierr = PetscSectionAddDof(network->DofSection,p,nvar);CHKERRQ(ierr); 1168 if (!compvalue) PetscFunctionReturn(0); 1169 1170 ierr = DMNetworkIsSharedVertex(dm,p,&sharedv);CHKERRQ(ierr); 1171 if (sharedv) { 1172 PetscBool ghost; 1173 ierr = DMNetworkIsGhostVertex(dm,p,&ghost);CHKERRQ(ierr); 1174 if (ghost) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Adding a component at a leaf(ghost) shared vertex is not supported"); 1175 } 1176 1177 if (compnum == MAX_DATA_AT_POINT) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components at a point exceeds the max %D",MAX_DATA_AT_POINT); 1178 1179 header->size[compnum] = component->size; 1180 ierr = PetscSectionAddDof(network->DataSection,p,component->size);CHKERRQ(ierr); 1181 header->key[compnum] = componentkey; 1182 if (compnum != 0) header->offset[compnum] = header->offset[compnum-1] + header->size[compnum-1]; 1183 cvalue->data[compnum] = (void*)compvalue; 1184 1185 /* variables */ 1186 header->nvar[compnum] += nvar; 1187 if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1]; 1188 1189 header->ndata++; 1190 PetscFunctionReturn(0); 1191 } 1192 1193 /*@ 1194 DMNetworkGetComponent - Gets the component key, the component data, and the number of variables at a given network point 1195 1196 Not Collective 1197 1198 Input Parameters: 1199 + dm - the DMNetwork object 1200 . p - vertex/edge point 1201 . compnum - component number; use ALL_COMPONENTS if sum up all the components 1202 1203 Output Parameters: 1204 + compkey - the key obtained when registering the component (use NULL if not required) 1205 . component - the component data (use NULL if not required) 1206 - nvar - number of variables (use NULL if not required) 1207 1208 Level: beginner 1209 1210 .seealso: DMNetworkAddComponent(), DMNetworkGetNumComponents() 1211 @*/ 1212 PetscErrorCode DMNetworkGetComponent(DM dm,PetscInt p,PetscInt compnum,PetscInt *compkey,void **component,PetscInt *nvar) 1213 { 1214 PetscErrorCode ierr; 1215 DM_Network *network = (DM_Network*)dm->data; 1216 PetscInt offset = 0; 1217 DMNetworkComponentHeader header; 1218 1219 PetscFunctionBegin; 1220 if (compnum == ALL_COMPONENTS) { 1221 ierr = PetscSectionGetDof(network->DofSection,p,nvar);CHKERRQ(ierr); 1222 PetscFunctionReturn(0); 1223 } 1224 1225 ierr = PetscSectionGetOffset(network->DataSection,p,&offset);CHKERRQ(ierr); 1226 header = (DMNetworkComponentHeader)(network->componentdataarray+offset);CHKERRQ(ierr); 1227 1228 if (compnum >= 0) { 1229 if (compkey) *compkey = header->key[compnum]; 1230 if (component) { 1231 offset += network->dataheadersize+header->offset[compnum]; 1232 *component = network->componentdataarray+offset; 1233 } 1234 } 1235 1236 if (nvar) *nvar = header->nvar[compnum]; 1237 PetscFunctionReturn(0); 1238 } 1239 1240 /* 1241 Sets up the array that holds the data for all components and its associated section. 1242 It copies the data for all components in a contiguous array called componentdataarray. The component data is stored pointwise with an additional header (metadata) stored for each point. The header has metadata information such as number of components at each point, number of variables for each component, offsets for the components data, etc. 1243 */ 1244 PetscErrorCode DMNetworkComponentSetUp(DM dm) 1245 { 1246 PetscErrorCode ierr; 1247 DM_Network *network = (DM_Network*)dm->data; 1248 PetscInt arr_size,p,offset,offsetp,ncomp,i; 1249 MPI_Comm comm; 1250 PetscMPIInt size,rank; 1251 DMNetworkComponentHeader header; 1252 DMNetworkComponentValue cvalue; 1253 DMNetworkComponentGenericDataType *componentdataarray; 1254 1255 PetscFunctionBegin; 1256 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1257 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 1258 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 1259 #if 0 1260 //------------- new 1261 if (size > 1) { /* Sync nvar at shared vertices for all processes */ 1262 PetscSF sf = network->plex->sf; 1263 const PetscInt *degree; 1264 PetscInt i,nleaves_total,*indata,*outdata,nroots,nleaves,nsv,p,ncomp; 1265 const PetscInt *svtx; 1266 PetscBool ghost; 1267 1268 ierr = PetscSFGetGraph(sf,&nroots,&nleaves,NULL,NULL);CHKERRQ(ierr); 1269 ierr = PetscSFComputeDegreeBegin(sf,°ree);CHKERRQ(ierr); 1270 ierr = PetscSFComputeDegreeEnd(sf,°ree);CHKERRQ(ierr); 1271 nleaves_total=0; 1272 for (i=0; i<nroots; i++) nleaves_total += degree[i]; 1273 printf("[%d] nleaves_total %d\n",rank,nleaves_total); 1274 MPI_Barrier(comm); 1275 1276 ierr = PetscCalloc2(nleaves_total,&indata,nleaves,&outdata);CHKERRQ(ierr); 1277 1278 /* Leaves copy user's ncomp to outdata */ 1279 ierr = DMNetworkGetSharedVertices(dm,&nsv,&svtx);CHKERRQ(ierr); 1280 for (i=0; i<nsv; i++) { 1281 p = svtx[i]; 1282 ierr = DMNetworkIsGhostVertex(dm,p,&ghost);CHKERRQ(ierr); 1283 if (!ghost) continue; 1284 1285 header = &network->header[p]; 1286 ncomp = header->ndata; 1287 printf("[%d] leaf has ncomp %d\n",rank,ncomp); 1288 outdata[p] = ncomp; 1289 } 1290 1291 /* Roots gather ncomp from leaves */ 1292 ierr = PetscSFGatherBegin(sf,MPIU_INT,outdata,indata);CHKERRQ(ierr); 1293 ierr = PetscSFGatherEnd(sf,MPIU_INT,outdata,indata);CHKERRQ(ierr); 1294 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"## Gathered data at multi-roots from leaves\n");CHKERRQ(ierr); 1295 ierr = PetscIntView(nleaves_total,indata,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1296 1297 ierr = PetscFree2(indata,outdata);CHKERRQ(ierr); 1298 } 1299 //---------------------- 1300 #endif 1301 1302 ierr = PetscSectionSetUp(network->DataSection);CHKERRQ(ierr); 1303 ierr = PetscSectionGetStorageSize(network->DataSection,&arr_size);CHKERRQ(ierr); 1304 ierr = PetscMalloc1(arr_size,&network->componentdataarray);CHKERRQ(ierr); 1305 componentdataarray = network->componentdataarray; 1306 for (p = network->pStart; p < network->pEnd; p++) { 1307 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 1308 /* Copy header */ 1309 header = &network->header[p]; 1310 ierr = PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));CHKERRQ(ierr); 1311 /* Copy data */ 1312 cvalue = &network->cvalue[p]; 1313 ncomp = header->ndata; 1314 1315 for (i = 0; i < ncomp; i++) { 1316 offset = offsetp + network->dataheadersize + header->offset[i]; 1317 ierr = PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));CHKERRQ(ierr); 1318 } 1319 } 1320 PetscFunctionReturn(0); 1321 } 1322 1323 /* Sets up the section for dofs. This routine is called during DMSetUp() */ 1324 static PetscErrorCode DMNetworkVariablesSetUp(DM dm) 1325 { 1326 PetscErrorCode ierr; 1327 DM_Network *network = (DM_Network*)dm->data; 1328 MPI_Comm comm; 1329 PetscMPIInt size; 1330 1331 PetscFunctionBegin; 1332 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1333 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 1334 1335 if (size > 1) { /* Sync nvar at shared vertices for all processes */ 1336 PetscSF sf = network->plex->sf; 1337 PetscInt *local_nvar, *remote_nvar,nroots,nleaves,p=-1,i,nsv; 1338 const PetscInt *svtx; 1339 PetscBool ghost; 1340 /* 1341 PetscMPIInt rank; 1342 const PetscInt *ilocal; 1343 const PetscSFNode *iremote; 1344 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 1345 ierr = PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 1346 */ 1347 ierr = PetscSFGetGraph(sf,&nroots,&nleaves,NULL,NULL);CHKERRQ(ierr); 1348 ierr = PetscCalloc2(nroots,&local_nvar,nroots,&remote_nvar);CHKERRQ(ierr); 1349 1350 /* Leaves copy user's nvar to local_nvar */ 1351 ierr = DMNetworkGetSharedVertices(dm,&nsv,&svtx);CHKERRQ(ierr); 1352 for (i=0; i<nsv; i++) { 1353 p = svtx[i]; 1354 ierr = DMNetworkIsGhostVertex(dm,p,&ghost);CHKERRQ(ierr); 1355 if (!ghost) continue; 1356 ierr = PetscSectionGetDof(network->DofSection,p,&local_nvar[p]);CHKERRQ(ierr); 1357 /* printf("[%d] Before SFReduce: leaf local_nvar[%d] = %d\n",rank,p,local_nvar[p]); */ 1358 } 1359 1360 /* Leaves add local_nvar to root remote_nvar */ 1361 ierr = PetscSFReduceBegin(sf, MPIU_INT, local_nvar, remote_nvar, MPI_SUM);CHKERRQ(ierr); 1362 ierr = PetscSFReduceEnd(sf, MPIU_INT, local_nvar, remote_nvar, MPI_SUM);CHKERRQ(ierr); 1363 /* printf("[%d] remote_nvar[%d] = %d\n",rank,p,remote_nvar[p]); */ 1364 1365 /* Update roots' local_nvar */ 1366 for (i=0; i<nsv; i++) { 1367 p = svtx[i]; 1368 ierr = DMNetworkIsGhostVertex(dm,p,&ghost);CHKERRQ(ierr); 1369 if (ghost) continue; 1370 1371 ierr = PetscSectionAddDof(network->DofSection,p,remote_nvar[p]);CHKERRQ(ierr); 1372 ierr = PetscSectionGetDof(network->DofSection,p,&local_nvar[p]);CHKERRQ(ierr); 1373 /* printf("[%d] After SFReduce: root local_nvar[%d] = %d\n",rank,p,local_nvar[p]); */ 1374 } 1375 1376 /* Roots Bcast nvar to leaves */ 1377 ierr = PetscSFBcastBegin(sf, MPIU_INT, local_nvar, remote_nvar);CHKERRQ(ierr); 1378 ierr = PetscSFBcastEnd(sf, MPIU_INT, local_nvar, remote_nvar);CHKERRQ(ierr); 1379 1380 /* Leaves reset receved/remote nvar to dm */ 1381 for (i=0; i<nsv; i++) { 1382 ierr = DMNetworkIsGhostVertex(dm,p,&ghost);CHKERRQ(ierr); 1383 if (!ghost) continue; 1384 p = svtx[i]; 1385 /* printf("[%d] leaf reset nvar %d at p= %d \n",rank,remote_nvar[p],p); */ 1386 /* DMNetworkSetNumVariables(dm,p,remote_nvar[p]) */ 1387 ierr = PetscSectionSetDof(network->DofSection,p,remote_nvar[p]);CHKERRQ(ierr); 1388 } 1389 1390 ierr = PetscFree2(local_nvar,remote_nvar);CHKERRQ(ierr); 1391 } 1392 1393 ierr = PetscSectionSetUp(network->DofSection);CHKERRQ(ierr); 1394 PetscFunctionReturn(0); 1395 } 1396 1397 /* Get a subsection from a range of points */ 1398 static PetscErrorCode DMNetworkGetSubSection_private(PetscSection master,PetscInt pstart,PetscInt pend,PetscSection *subsection) 1399 { 1400 PetscErrorCode ierr; 1401 PetscInt i, nvar; 1402 1403 PetscFunctionBegin; 1404 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);CHKERRQ(ierr); 1405 ierr = PetscSectionSetChart(*subsection, 0, pend - pstart);CHKERRQ(ierr); 1406 for (i = pstart; i < pend; i++) { 1407 ierr = PetscSectionGetDof(master,i,&nvar);CHKERRQ(ierr); 1408 ierr = PetscSectionSetDof(*subsection, i - pstart, nvar);CHKERRQ(ierr); 1409 } 1410 1411 ierr = PetscSectionSetUp(*subsection);CHKERRQ(ierr); 1412 PetscFunctionReturn(0); 1413 } 1414 1415 /* Create a submap of points with a GlobalToLocal structure */ 1416 static PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map) 1417 { 1418 PetscErrorCode ierr; 1419 PetscInt i, *subpoints; 1420 1421 PetscFunctionBegin; 1422 /* Create index sets to map from "points" to "subpoints" */ 1423 ierr = PetscMalloc1(pend - pstart, &subpoints);CHKERRQ(ierr); 1424 for (i = pstart; i < pend; i++) { 1425 subpoints[i - pstart] = i; 1426 } 1427 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);CHKERRQ(ierr); 1428 ierr = PetscFree(subpoints);CHKERRQ(ierr); 1429 PetscFunctionReturn(0); 1430 } 1431 1432 /*@ 1433 DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute 1434 1435 Collective on dm 1436 1437 Input Parameters: 1438 . dm - the DMNetworkObject 1439 1440 Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are: 1441 1442 points = [0 1 2 3 4 5 6] 1443 1444 where edges = [0,1,2,3] and vertices = [4,5,6]. The new orderings will be specific to the subset (i.e vertices = [0,1,2] <- [4,5,6]). 1445 1446 With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset. 1447 1448 Level: intermediate 1449 1450 @*/ 1451 PetscErrorCode DMNetworkAssembleGraphStructures(DM dm) 1452 { 1453 PetscErrorCode ierr; 1454 MPI_Comm comm; 1455 PetscMPIInt rank, size; 1456 DM_Network *network = (DM_Network*)dm->data; 1457 1458 PetscFunctionBegin; 1459 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1460 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1461 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1462 1463 /* Create maps for vertices and edges */ 1464 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 1465 ierr = DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);CHKERRQ(ierr); 1466 1467 /* Create local sub-sections */ 1468 ierr = DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);CHKERRQ(ierr); 1469 ierr = DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);CHKERRQ(ierr); 1470 1471 if (size > 1) { 1472 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 1473 1474 ierr = PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);CHKERRQ(ierr); 1475 ierr = PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);CHKERRQ(ierr); 1476 ierr = PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);CHKERRQ(ierr); 1477 } else { 1478 /* create structures for vertex */ 1479 ierr = PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);CHKERRQ(ierr); 1480 /* create structures for edge */ 1481 ierr = PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);CHKERRQ(ierr); 1482 } 1483 1484 /* Add viewers */ 1485 ierr = PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");CHKERRQ(ierr); 1486 ierr = PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");CHKERRQ(ierr); 1487 ierr = PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");CHKERRQ(ierr); 1488 ierr = PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");CHKERRQ(ierr); 1489 PetscFunctionReturn(0); 1490 } 1491 1492 /*@ 1493 DMNetworkDistribute - Distributes the network and moves associated component data 1494 1495 Collective 1496 1497 Input Parameter: 1498 + DM - the DMNetwork object 1499 - overlap - the overlap of partitions, 0 is the default 1500 1501 Notes: 1502 Distributes the network with <overlap>-overlapping partitioning of the edges. 1503 1504 Level: intermediate 1505 1506 .seealso: DMNetworkCreate() 1507 @*/ 1508 PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap) 1509 { 1510 MPI_Comm comm; 1511 PetscErrorCode ierr; 1512 PetscMPIInt size; 1513 DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data); 1514 DM_Network *newDMnetwork; 1515 PetscSF pointsf=NULL; 1516 DM newDM; 1517 PetscInt j,e,v,offset,*subnetvtx,nsubnet,gidx,svtx_idx,nv; 1518 PetscInt to_net,from_net,*svto; 1519 PetscPartitioner part; 1520 DMNetworkComponentHeader header; 1521 1522 PetscFunctionBegin; 1523 ierr = PetscObjectGetComm((PetscObject)*dm,&comm);CHKERRQ(ierr); 1524 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1525 if (size == 1) PetscFunctionReturn(0); 1526 1527 /* This routine moves the component data to the appropriate processors. It makes use of the DataSection and the componentdataarray to move the component data to appropriate processors and returns a new DataSection and new componentdataarray. */ 1528 ierr = DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);CHKERRQ(ierr); 1529 newDMnetwork = (DM_Network*)newDM->data; 1530 newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType); 1531 1532 /* Enable runtime options for petscpartitioner */ 1533 ierr = DMPlexGetPartitioner(oldDMnetwork->plex,&part);CHKERRQ(ierr); 1534 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 1535 1536 /* Distribute plex dm */ 1537 ierr = DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);CHKERRQ(ierr); 1538 1539 /* Distribute dof section */ 1540 ierr = PetscSectionCreate(comm,&newDMnetwork->DofSection);CHKERRQ(ierr); 1541 ierr = PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);CHKERRQ(ierr); 1542 1543 /* Distribute data and associated section */ 1544 ierr = PetscSectionCreate(comm,&newDMnetwork->DataSection);CHKERRQ(ierr); 1545 ierr = DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);CHKERRQ(ierr); 1546 1547 ierr = PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);CHKERRQ(ierr); 1548 ierr = DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);CHKERRQ(ierr); 1549 ierr = DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);CHKERRQ(ierr); 1550 newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart; 1551 newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart; 1552 newDMnetwork->NVertices = oldDMnetwork->NVertices; 1553 newDMnetwork->NEdges = oldDMnetwork->NEdges; 1554 newDMnetwork->svtable = oldDMnetwork->svtable; 1555 oldDMnetwork->svtable = NULL; 1556 1557 /* Set Dof section as the section for dm */ 1558 ierr = DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);CHKERRQ(ierr); 1559 ierr = DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);CHKERRQ(ierr); 1560 1561 /* Set up subnetwork info in the newDM */ 1562 newDMnetwork->Nsubnet = oldDMnetwork->Nsubnet; 1563 newDMnetwork->Nsvtx = oldDMnetwork->Nsvtx; 1564 oldDMnetwork->Nsvtx = 0; 1565 newDMnetwork->svtx = oldDMnetwork->svtx; /* global vertices! */ 1566 oldDMnetwork->svtx = NULL; 1567 ierr = PetscCalloc1(newDMnetwork->Nsubnet,&newDMnetwork->subnet);CHKERRQ(ierr); 1568 1569 /* Copy over the global number of vertices and edges in each subnetwork. 1570 Note: these are calculated in DMNetworkLayoutSetUp() 1571 */ 1572 nsubnet = newDMnetwork->Nsubnet; 1573 for (j = 0; j < nsubnet; j++) { 1574 newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx; 1575 newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge; 1576 } 1577 1578 /* Get local nedges and nvtx for subnetworks */ 1579 for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) { 1580 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);CHKERRQ(ierr); 1581 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1582 newDMnetwork->subnet[header->subnetid].nedge++; 1583 } 1584 1585 for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) { 1586 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);CHKERRQ(ierr); 1587 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1588 1589 /* shared vertices: use gidx = header->index to check if v is a shared vertex */ 1590 gidx = header->index; 1591 ierr = PetscTableFind(newDMnetwork->svtable,gidx+1,&svtx_idx);CHKERRQ(ierr); 1592 svtx_idx--; 1593 1594 if (svtx_idx < 0) { /* not a shared vertex */ 1595 newDMnetwork->subnet[header->subnetid].nvtx++; 1596 } else { /* a shared vertex belongs to more than one subnetworks, it is being counted by multiple subnets */ 1597 from_net = newDMnetwork->svtx[svtx_idx].sv[0]; 1598 newDMnetwork->subnet[from_net].nvtx++; 1599 for (j=1; j<newDMnetwork->svtx[svtx_idx].n; j++) { 1600 svto = newDMnetwork->svtx[svtx_idx].sv + 2*j; 1601 to_net = svto[0]; 1602 newDMnetwork->subnet[to_net].nvtx++; 1603 } 1604 } 1605 } 1606 1607 /* Get total local nvtx for subnetworks */ 1608 nv = 0; 1609 for (j=0; j<nsubnet; j++) nv += newDMnetwork->subnet[j].nvtx; 1610 nv += newDMnetwork->Nsvtx; 1611 1612 /* Now create the vertices and edge arrays for the subnetworks */ 1613 ierr = PetscCalloc1(nv,&subnetvtx);CHKERRQ(ierr); 1614 newDMnetwork->subnetvtx = subnetvtx; 1615 1616 for (j=0; j<nsubnet; j++) { 1617 ierr = PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);CHKERRQ(ierr); 1618 newDMnetwork->subnet[j].vertices = subnetvtx; 1619 subnetvtx += newDMnetwork->subnet[j].nvtx; 1620 1621 /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop. These get updated when the vertices and edges are added. */ 1622 newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0; 1623 } 1624 newDMnetwork->svertices = subnetvtx; 1625 1626 /* Set the edges and vertices in each subnetwork */ 1627 for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) { 1628 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);CHKERRQ(ierr); 1629 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1630 newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e; 1631 } 1632 1633 nv = 0; 1634 for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) { 1635 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);CHKERRQ(ierr); 1636 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1637 1638 /* coupling vertices: use gidx = header->index to check if v is a coupling vertex */ 1639 ierr = PetscTableFind(newDMnetwork->svtable,header->index+1,&svtx_idx);CHKERRQ(ierr); 1640 svtx_idx--; 1641 if (svtx_idx < 0) { 1642 newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v; 1643 } else { /* a shared vertex */ 1644 newDMnetwork->svertices[nv++] = v; 1645 1646 from_net = newDMnetwork->svtx[svtx_idx].sv[0]; 1647 newDMnetwork->subnet[from_net].vertices[newDMnetwork->subnet[from_net].nvtx++] = v; 1648 1649 for (j=1; j<newDMnetwork->svtx[svtx_idx].n; j++) { 1650 svto = newDMnetwork->svtx[svtx_idx].sv + 2*j; 1651 to_net = svto[0]; 1652 newDMnetwork->subnet[to_net].vertices[newDMnetwork->subnet[to_net].nvtx++] = v; 1653 } 1654 } 1655 } 1656 newDMnetwork->nsvtx = nv; /* num of local shared vertices */ 1657 1658 newDM->setupcalled = (*dm)->setupcalled; 1659 newDMnetwork->distributecalled = PETSC_TRUE; 1660 1661 /* Free spaces */ 1662 ierr = PetscSFDestroy(&pointsf);CHKERRQ(ierr); 1663 ierr = DMDestroy(dm);CHKERRQ(ierr); 1664 1665 *dm = newDM; 1666 PetscFunctionReturn(0); 1667 } 1668 1669 /*@C 1670 PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering 1671 1672 Collective 1673 1674 Input Parameters: 1675 + masterSF - the original SF structure 1676 - map - a ISLocalToGlobal mapping that contains the subset of points 1677 1678 Output Parameters: 1679 . subSF - a subset of the masterSF for the desired subset. 1680 1681 Level: intermediate 1682 @*/ 1683 PetscErrorCode PetscSFGetSubSF(PetscSF mastersf,ISLocalToGlobalMapping map,PetscSF *subSF) 1684 { 1685 PetscErrorCode ierr; 1686 PetscInt nroots, nleaves, *ilocal_sub; 1687 PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0; 1688 PetscInt *local_points, *remote_points; 1689 PetscSFNode *iremote_sub; 1690 const PetscInt *ilocal; 1691 const PetscSFNode *iremote; 1692 1693 PetscFunctionBegin; 1694 ierr = PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 1695 1696 /* Look for leaves that pertain to the subset of points. Get the local ordering */ 1697 ierr = PetscMalloc1(nleaves,&ilocal_map);CHKERRQ(ierr); 1698 ierr = ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);CHKERRQ(ierr); 1699 for (i = 0; i < nleaves; i++) { 1700 if (ilocal_map[i] != -1) nleaves_sub += 1; 1701 } 1702 /* Re-number ilocal with subset numbering. Need information from roots */ 1703 ierr = PetscMalloc2(nroots,&local_points,nroots,&remote_points);CHKERRQ(ierr); 1704 for (i = 0; i < nroots; i++) local_points[i] = i; 1705 ierr = ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);CHKERRQ(ierr); 1706 ierr = PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);CHKERRQ(ierr); 1707 ierr = PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);CHKERRQ(ierr); 1708 /* Fill up graph using local (that is, local to the subset) numbering. */ 1709 ierr = PetscMalloc1(nleaves_sub,&ilocal_sub);CHKERRQ(ierr); 1710 ierr = PetscMalloc1(nleaves_sub,&iremote_sub);CHKERRQ(ierr); 1711 nleaves_sub = 0; 1712 for (i = 0; i < nleaves; i++) { 1713 if (ilocal_map[i] != -1) { 1714 ilocal_sub[nleaves_sub] = ilocal_map[i]; 1715 iremote_sub[nleaves_sub].rank = iremote[i].rank; 1716 iremote_sub[nleaves_sub].index = remote_points[ilocal[i]]; 1717 nleaves_sub += 1; 1718 } 1719 } 1720 ierr = PetscFree2(local_points,remote_points);CHKERRQ(ierr); 1721 ierr = ISLocalToGlobalMappingGetSize(map,&nroots_sub);CHKERRQ(ierr); 1722 1723 /* Create new subSF */ 1724 ierr = PetscSFCreate(PETSC_COMM_WORLD,subSF);CHKERRQ(ierr); 1725 ierr = PetscSFSetFromOptions(*subSF);CHKERRQ(ierr); 1726 ierr = PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);CHKERRQ(ierr); 1727 ierr = PetscFree(ilocal_map);CHKERRQ(ierr); 1728 ierr = PetscFree(iremote_sub);CHKERRQ(ierr); 1729 PetscFunctionReturn(0); 1730 } 1731 1732 /*@C 1733 DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point 1734 1735 Not Collective 1736 1737 Input Parameters: 1738 + dm - the DMNetwork object 1739 - p - the vertex point 1740 1741 Output Parameters: 1742 + nedges - number of edges connected to this vertex point 1743 - edges - list of edge points 1744 1745 Level: beginner 1746 1747 Fortran Notes: 1748 Since it returns an array, this routine is only available in Fortran 90, and you must 1749 include petsc.h90 in your code. 1750 1751 .seealso: DMNetworkCreate(), DMNetworkGetConnectedVertices() 1752 @*/ 1753 PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[]) 1754 { 1755 PetscErrorCode ierr; 1756 DM_Network *network = (DM_Network*)dm->data; 1757 1758 PetscFunctionBegin; 1759 ierr = DMPlexGetSupportSize(network->plex,vertex,nedges);CHKERRQ(ierr); 1760 ierr = DMPlexGetSupport(network->plex,vertex,edges);CHKERRQ(ierr); 1761 PetscFunctionReturn(0); 1762 } 1763 1764 /*@C 1765 DMNetworkGetConnectedVertices - Return the connected vertices for this edge point 1766 1767 Not Collective 1768 1769 Input Parameters: 1770 + dm - the DMNetwork object 1771 - p - the edge point 1772 1773 Output Parameters: 1774 . vertices - vertices connected to this edge 1775 1776 Level: beginner 1777 1778 Fortran Notes: 1779 Since it returns an array, this routine is only available in Fortran 90, and you must 1780 include petsc.h90 in your code. 1781 1782 .seealso: DMNetworkCreate(), DMNetworkGetSupportingEdges() 1783 @*/ 1784 PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[]) 1785 { 1786 PetscErrorCode ierr; 1787 DM_Network *network = (DM_Network*)dm->data; 1788 1789 PetscFunctionBegin; 1790 ierr = DMPlexGetCone(network->plex,edge,vertices);CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /*@ 1795 DMNetworkIsSharedVertex - Returns TRUE if the vertex is shared by subnetworks 1796 1797 Not Collective 1798 1799 Input Parameters: 1800 + dm - the DMNetwork object 1801 - p - the vertex point 1802 1803 Output Parameter: 1804 . flag - TRUE if the vertex is shared by subnetworks 1805 1806 Level: beginner 1807 1808 .seealso: DMNetworkAddSharedVertices(), DMNetworkIsGhostVertex() 1809 @*/ 1810 PetscErrorCode DMNetworkIsSharedVertex(DM dm,PetscInt p,PetscBool *flag) 1811 { 1812 PetscErrorCode ierr; 1813 PetscInt i; 1814 1815 PetscFunctionBegin; 1816 *flag = PETSC_FALSE; 1817 1818 if (dm->setupcalled) { /* DMNetworkGetGlobalVertexIndex() requires DMSetUp() be called */ 1819 DM_Network *network = (DM_Network*)dm->data; 1820 PetscInt gidx; 1821 ierr = DMNetworkGetGlobalVertexIndex(dm,p,&gidx);CHKERRQ(ierr); 1822 ierr = PetscTableFind(network->svtable,gidx+1,&i);CHKERRQ(ierr); 1823 if (i) *flag = PETSC_TRUE; 1824 } else { /* would be removed? */ 1825 PetscInt nv; 1826 const PetscInt *vtx; 1827 ierr = DMNetworkGetSharedVertices(dm,&nv,&vtx);CHKERRQ(ierr); 1828 for (i=0; i<nv; i++) { 1829 if (p == vtx[i]) { 1830 *flag = PETSC_TRUE; 1831 break; 1832 } 1833 } 1834 } 1835 PetscFunctionReturn(0); 1836 } 1837 1838 /*@ 1839 DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex 1840 1841 Not Collective 1842 1843 Input Parameters: 1844 + dm - the DMNetwork object 1845 - p - the vertex point 1846 1847 Output Parameter: 1848 . isghost - TRUE if the vertex is a ghost point 1849 1850 Level: beginner 1851 1852 .seealso: DMNetworkGetConnectedVertices(), DMNetworkGetVertexRange(), DMNetworkIsSharedVertex() 1853 @*/ 1854 PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost) 1855 { 1856 PetscErrorCode ierr; 1857 DM_Network *network = (DM_Network*)dm->data; 1858 PetscInt offsetg; 1859 PetscSection sectiong; 1860 1861 PetscFunctionBegin; 1862 *isghost = PETSC_FALSE; 1863 ierr = DMGetGlobalSection(network->plex,§iong);CHKERRQ(ierr); 1864 ierr = PetscSectionGetOffset(sectiong,p,&offsetg);CHKERRQ(ierr); 1865 if (offsetg < 0) *isghost = PETSC_TRUE; 1866 PetscFunctionReturn(0); 1867 } 1868 1869 PetscErrorCode DMSetUp_Network(DM dm) 1870 { 1871 PetscErrorCode ierr; 1872 DM_Network *network=(DM_Network*)dm->data; 1873 1874 PetscFunctionBegin; 1875 ierr = DMNetworkComponentSetUp(dm);CHKERRQ(ierr); 1876 ierr = DMNetworkVariablesSetUp(dm);CHKERRQ(ierr); 1877 1878 ierr = DMSetLocalSection(network->plex,network->DofSection);CHKERRQ(ierr); 1879 ierr = DMGetGlobalSection(network->plex,&network->GlobalDofSection);CHKERRQ(ierr); 1880 1881 dm->setupcalled = PETSC_TRUE; 1882 ierr = DMViewFromOptions(dm,NULL,"-dm_view");CHKERRQ(ierr); 1883 PetscFunctionReturn(0); 1884 } 1885 1886 /*@ 1887 DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices 1888 -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)? 1889 1890 Collective 1891 1892 Input Parameters: 1893 + dm - the DMNetwork object 1894 . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges 1895 - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices 1896 1897 Level: intermediate 1898 1899 @*/ 1900 PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg) 1901 { 1902 DM_Network *network=(DM_Network*)dm->data; 1903 PetscErrorCode ierr; 1904 PetscInt nVertices = network->nVertices; 1905 1906 PetscFunctionBegin; 1907 network->userEdgeJacobian = eflg; 1908 network->userVertexJacobian = vflg; 1909 1910 if (eflg && !network->Je) { 1911 ierr = PetscCalloc1(3*network->nEdges,&network->Je);CHKERRQ(ierr); 1912 } 1913 1914 if (vflg && !network->Jv && nVertices) { 1915 PetscInt i,*vptr,nedges,vStart=network->vStart; 1916 PetscInt nedges_total; 1917 const PetscInt *edges; 1918 1919 /* count nvertex_total */ 1920 nedges_total = 0; 1921 ierr = PetscMalloc1(nVertices+1,&vptr);CHKERRQ(ierr); 1922 1923 vptr[0] = 0; 1924 for (i=0; i<nVertices; i++) { 1925 ierr = DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);CHKERRQ(ierr); 1926 nedges_total += nedges; 1927 vptr[i+1] = vptr[i] + 2*nedges + 1; 1928 } 1929 1930 ierr = PetscCalloc1(2*nedges_total+nVertices,&network->Jv);CHKERRQ(ierr); 1931 network->Jvptr = vptr; 1932 } 1933 PetscFunctionReturn(0); 1934 } 1935 1936 /*@ 1937 DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network 1938 1939 Not Collective 1940 1941 Input Parameters: 1942 + dm - the DMNetwork object 1943 . p - the edge point 1944 - J - array (size = 3) of Jacobian submatrices for this edge point: 1945 J[0]: this edge 1946 J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices() 1947 1948 Level: advanced 1949 1950 .seealso: DMNetworkVertexSetMatrix() 1951 @*/ 1952 PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[]) 1953 { 1954 DM_Network *network=(DM_Network*)dm->data; 1955 1956 PetscFunctionBegin; 1957 if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix"); 1958 1959 if (J) { 1960 network->Je[3*p] = J[0]; 1961 network->Je[3*p+1] = J[1]; 1962 network->Je[3*p+2] = J[2]; 1963 } 1964 PetscFunctionReturn(0); 1965 } 1966 1967 /*@ 1968 DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network 1969 1970 Not Collective 1971 1972 Input Parameters: 1973 + dm - The DMNetwork object 1974 . p - the vertex point 1975 - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point: 1976 J[0]: this vertex 1977 J[1+2*i]: i-th supporting edge 1978 J[1+2*i+1]: i-th connected vertex 1979 1980 Level: advanced 1981 1982 .seealso: DMNetworkEdgeSetMatrix() 1983 @*/ 1984 PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[]) 1985 { 1986 PetscErrorCode ierr; 1987 DM_Network *network=(DM_Network*)dm->data; 1988 PetscInt i,*vptr,nedges,vStart=network->vStart; 1989 const PetscInt *edges; 1990 1991 PetscFunctionBegin; 1992 if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix"); 1993 1994 if (J) { 1995 vptr = network->Jvptr; 1996 network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */ 1997 1998 /* Set Jacobian for each supporting edge and connected vertex */ 1999 ierr = DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);CHKERRQ(ierr); 2000 for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i]; 2001 } 2002 PetscFunctionReturn(0); 2003 } 2004 2005 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 2006 { 2007 PetscErrorCode ierr; 2008 PetscInt j; 2009 PetscScalar val=(PetscScalar)ncols; 2010 2011 PetscFunctionBegin; 2012 if (!ghost) { 2013 for (j=0; j<nrows; j++) { 2014 ierr = VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 2015 } 2016 } else { 2017 for (j=0; j<nrows; j++) { 2018 ierr = VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 2019 } 2020 } 2021 PetscFunctionReturn(0); 2022 } 2023 2024 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 2025 { 2026 PetscErrorCode ierr; 2027 PetscInt j,ncols_u; 2028 PetscScalar val; 2029 2030 PetscFunctionBegin; 2031 if (!ghost) { 2032 for (j=0; j<nrows; j++) { 2033 ierr = MatGetRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 2034 val = (PetscScalar)ncols_u; 2035 ierr = VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 2036 ierr = MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 2037 } 2038 } else { 2039 for (j=0; j<nrows; j++) { 2040 ierr = MatGetRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 2041 val = (PetscScalar)ncols_u; 2042 ierr = VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 2043 ierr = MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 2044 } 2045 } 2046 PetscFunctionReturn(0); 2047 } 2048 2049 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 2050 { 2051 PetscErrorCode ierr; 2052 2053 PetscFunctionBegin; 2054 if (Ju) { 2055 ierr = MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);CHKERRQ(ierr); 2056 } else { 2057 ierr = MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);CHKERRQ(ierr); 2058 } 2059 PetscFunctionReturn(0); 2060 } 2061 2062 PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 2063 { 2064 PetscErrorCode ierr; 2065 PetscInt j,*cols; 2066 PetscScalar *zeros; 2067 2068 PetscFunctionBegin; 2069 ierr = PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);CHKERRQ(ierr); 2070 for (j=0; j<ncols; j++) cols[j] = j+ cstart; 2071 ierr = MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);CHKERRQ(ierr); 2072 ierr = PetscFree2(cols,zeros);CHKERRQ(ierr); 2073 PetscFunctionReturn(0); 2074 } 2075 2076 PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 2077 { 2078 PetscErrorCode ierr; 2079 PetscInt j,M,N,row,col,ncols_u; 2080 const PetscInt *cols; 2081 PetscScalar zero=0.0; 2082 2083 PetscFunctionBegin; 2084 ierr = MatGetSize(Ju,&M,&N);CHKERRQ(ierr); 2085 if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N); 2086 2087 for (row=0; row<nrows; row++) { 2088 ierr = MatGetRow(Ju,row,&ncols_u,&cols,NULL);CHKERRQ(ierr); 2089 for (j=0; j<ncols_u; j++) { 2090 col = cols[j] + cstart; 2091 ierr = MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);CHKERRQ(ierr); 2092 } 2093 ierr = MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);CHKERRQ(ierr); 2094 } 2095 PetscFunctionReturn(0); 2096 } 2097 2098 PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 2099 { 2100 PetscErrorCode ierr; 2101 2102 PetscFunctionBegin; 2103 if (Ju) { 2104 ierr = MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 2105 } else { 2106 ierr = MatSetDenseblock_private(nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 2107 } 2108 PetscFunctionReturn(0); 2109 } 2110 2111 /* Creates a GlobalToLocal mapping with a Local and Global section. This is akin to the routine DMGetLocalToGlobalMapping but without the need of providing a dm. 2112 */ 2113 PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog) 2114 { 2115 PetscErrorCode ierr; 2116 PetscInt i,size,dof; 2117 PetscInt *glob2loc; 2118 2119 PetscFunctionBegin; 2120 ierr = PetscSectionGetStorageSize(localsec,&size);CHKERRQ(ierr); 2121 ierr = PetscMalloc1(size,&glob2loc);CHKERRQ(ierr); 2122 2123 for (i = 0; i < size; i++) { 2124 ierr = PetscSectionGetOffset(globalsec,i,&dof);CHKERRQ(ierr); 2125 dof = (dof >= 0) ? dof : -(dof + 1); 2126 glob2loc[i] = dof; 2127 } 2128 2129 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 2130 #if 0 2131 ierr = PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 2132 #endif 2133 PetscFunctionReturn(0); 2134 } 2135 2136 #include <petsc/private/matimpl.h> 2137 2138 PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J) 2139 { 2140 PetscErrorCode ierr; 2141 DM_Network *network = (DM_Network*)dm->data; 2142 PetscMPIInt rank, size; 2143 PetscInt eDof,vDof; 2144 Mat j11,j12,j21,j22,bA[2][2]; 2145 MPI_Comm comm; 2146 ISLocalToGlobalMapping eISMap,vISMap; 2147 2148 PetscFunctionBegin; 2149 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2150 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 2151 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 2152 2153 ierr = PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);CHKERRQ(ierr); 2154 ierr = PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);CHKERRQ(ierr); 2155 2156 ierr = MatCreate(comm, &j11);CHKERRQ(ierr); 2157 ierr = MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 2158 ierr = MatSetType(j11, MATMPIAIJ);CHKERRQ(ierr); 2159 2160 ierr = MatCreate(comm, &j12);CHKERRQ(ierr); 2161 ierr = MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);CHKERRQ(ierr); 2162 ierr = MatSetType(j12, MATMPIAIJ);CHKERRQ(ierr); 2163 2164 ierr = MatCreate(comm, &j21);CHKERRQ(ierr); 2165 ierr = MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 2166 ierr = MatSetType(j21, MATMPIAIJ);CHKERRQ(ierr); 2167 2168 ierr = MatCreate(comm, &j22);CHKERRQ(ierr); 2169 ierr = MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 2170 ierr = MatSetType(j22, MATMPIAIJ);CHKERRQ(ierr); 2171 2172 bA[0][0] = j11; 2173 bA[0][1] = j12; 2174 bA[1][0] = j21; 2175 bA[1][1] = j22; 2176 2177 ierr = CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);CHKERRQ(ierr); 2178 ierr = CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);CHKERRQ(ierr); 2179 2180 ierr = MatSetLocalToGlobalMapping(j11,eISMap,eISMap);CHKERRQ(ierr); 2181 ierr = MatSetLocalToGlobalMapping(j12,eISMap,vISMap);CHKERRQ(ierr); 2182 ierr = MatSetLocalToGlobalMapping(j21,vISMap,eISMap);CHKERRQ(ierr); 2183 ierr = MatSetLocalToGlobalMapping(j22,vISMap,vISMap);CHKERRQ(ierr); 2184 2185 ierr = MatSetUp(j11);CHKERRQ(ierr); 2186 ierr = MatSetUp(j12);CHKERRQ(ierr); 2187 ierr = MatSetUp(j21);CHKERRQ(ierr); 2188 ierr = MatSetUp(j22);CHKERRQ(ierr); 2189 2190 ierr = MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);CHKERRQ(ierr); 2191 ierr = MatSetUp(*J);CHKERRQ(ierr); 2192 ierr = MatNestSetVecType(*J,VECNEST);CHKERRQ(ierr); 2193 ierr = MatDestroy(&j11);CHKERRQ(ierr); 2194 ierr = MatDestroy(&j12);CHKERRQ(ierr); 2195 ierr = MatDestroy(&j21);CHKERRQ(ierr); 2196 ierr = MatDestroy(&j22);CHKERRQ(ierr); 2197 2198 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2199 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2200 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2201 2202 /* Free structures */ 2203 ierr = ISLocalToGlobalMappingDestroy(&eISMap);CHKERRQ(ierr); 2204 ierr = ISLocalToGlobalMappingDestroy(&vISMap);CHKERRQ(ierr); 2205 PetscFunctionReturn(0); 2206 } 2207 2208 PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J) 2209 { 2210 PetscErrorCode ierr; 2211 DM_Network *network = (DM_Network*)dm->data; 2212 PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize; 2213 PetscInt cstart,ncols,j,e,v; 2214 PetscBool ghost,ghost_vc,ghost2,isNest; 2215 Mat Juser; 2216 PetscSection sectionGlobal; 2217 PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */ 2218 const PetscInt *edges,*cone; 2219 MPI_Comm comm; 2220 MatType mtype; 2221 Vec vd_nz,vo_nz; 2222 PetscInt *dnnz,*onnz; 2223 PetscScalar *vdnz,*vonz; 2224 2225 PetscFunctionBegin; 2226 mtype = dm->mattype; 2227 ierr = PetscStrcmp(mtype,MATNEST,&isNest);CHKERRQ(ierr); 2228 if (isNest) { 2229 ierr = DMCreateMatrix_Network_Nest(dm,J);CHKERRQ(ierr); 2230 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 2231 PetscFunctionReturn(0); 2232 } 2233 2234 if (!network->userEdgeJacobian && !network->userVertexJacobian) { 2235 /* user does not provide Jacobian blocks */ 2236 ierr = DMCreateMatrix_Plex(network->plex,J);CHKERRQ(ierr); 2237 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 2238 PetscFunctionReturn(0); 2239 } 2240 2241 ierr = MatCreate(PetscObjectComm((PetscObject)dm),J);CHKERRQ(ierr); 2242 ierr = DMGetGlobalSection(network->plex,§ionGlobal);CHKERRQ(ierr); 2243 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);CHKERRQ(ierr); 2244 ierr = MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2245 2246 ierr = MatSetType(*J,MATAIJ);CHKERRQ(ierr); 2247 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 2248 2249 /* (1) Set matrix preallocation */ 2250 /*------------------------------*/ 2251 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2252 ierr = VecCreate(comm,&vd_nz);CHKERRQ(ierr); 2253 ierr = VecSetSizes(vd_nz,localSize,PETSC_DECIDE);CHKERRQ(ierr); 2254 ierr = VecSetFromOptions(vd_nz);CHKERRQ(ierr); 2255 ierr = VecSet(vd_nz,0.0);CHKERRQ(ierr); 2256 ierr = VecDuplicate(vd_nz,&vo_nz);CHKERRQ(ierr); 2257 2258 /* Set preallocation for edges */ 2259 /*-----------------------------*/ 2260 ierr = DMNetworkGetEdgeRange(dm,&eStart,&eEnd);CHKERRQ(ierr); 2261 2262 ierr = PetscMalloc1(localSize,&rows);CHKERRQ(ierr); 2263 for (e=eStart; e<eEnd; e++) { 2264 /* Get row indices */ 2265 ierr = DMNetworkGetGlobalVecOffset(dm,e,ALL_COMPONENTS,&rstart);CHKERRQ(ierr); 2266 ierr = PetscSectionGetDof(network->DofSection,e,&nrows);CHKERRQ(ierr); 2267 if (nrows) { 2268 for (j=0; j<nrows; j++) rows[j] = j + rstart; 2269 2270 /* Set preallocation for conntected vertices */ 2271 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 2272 for (v=0; v<2; v++) { 2273 ierr = PetscSectionGetDof(network->DofSection,cone[v],&ncols);CHKERRQ(ierr); 2274 2275 if (network->Je) { 2276 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 2277 } else Juser = NULL; 2278 ierr = DMNetworkIsGhostVertex(dm,cone[v],&ghost);CHKERRQ(ierr); 2279 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 2280 } 2281 2282 /* Set preallocation for edge self */ 2283 cstart = rstart; 2284 if (network->Je) { 2285 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 2286 } else Juser = NULL; 2287 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 2288 } 2289 } 2290 2291 /* Set preallocation for vertices */ 2292 /*--------------------------------*/ 2293 ierr = DMNetworkGetVertexRange(dm,&vStart,&vEnd);CHKERRQ(ierr); 2294 if (vEnd - vStart) vptr = network->Jvptr; 2295 2296 for (v=vStart; v<vEnd; v++) { 2297 /* Get row indices */ 2298 ierr = DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&rstart);CHKERRQ(ierr); 2299 ierr = PetscSectionGetDof(network->DofSection,v,&nrows);CHKERRQ(ierr); 2300 if (!nrows) continue; 2301 2302 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 2303 if (ghost) { 2304 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 2305 } else { 2306 rows_v = rows; 2307 } 2308 2309 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 2310 2311 /* Get supporting edges and connected vertices */ 2312 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 2313 2314 for (e=0; e<nedges; e++) { 2315 /* Supporting edges */ 2316 ierr = DMNetworkGetGlobalVecOffset(dm,edges[e],ALL_COMPONENTS,&cstart);CHKERRQ(ierr); 2317 ierr = PetscSectionGetDof(network->DofSection,edges[e],&ncols);CHKERRQ(ierr); 2318 2319 if (network->Jv) { 2320 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 2321 } else Juser = NULL; 2322 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 2323 2324 /* Connected vertices */ 2325 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 2326 vc = (v == cone[0]) ? cone[1]:cone[0]; 2327 ierr = DMNetworkIsGhostVertex(dm,vc,&ghost_vc);CHKERRQ(ierr); 2328 2329 ierr = PetscSectionGetDof(network->DofSection,vc,&ncols);CHKERRQ(ierr); 2330 2331 if (network->Jv) { 2332 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 2333 } else Juser = NULL; 2334 if (ghost_vc||ghost) { 2335 ghost2 = PETSC_TRUE; 2336 } else { 2337 ghost2 = PETSC_FALSE; 2338 } 2339 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);CHKERRQ(ierr); 2340 } 2341 2342 /* Set preallocation for vertex self */ 2343 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 2344 if (!ghost) { 2345 ierr = DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&cstart);CHKERRQ(ierr); 2346 if (network->Jv) { 2347 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 2348 } else Juser = NULL; 2349 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 2350 } 2351 if (ghost) { 2352 ierr = PetscFree(rows_v);CHKERRQ(ierr); 2353 } 2354 } 2355 2356 ierr = VecAssemblyBegin(vd_nz);CHKERRQ(ierr); 2357 ierr = VecAssemblyBegin(vo_nz);CHKERRQ(ierr); 2358 2359 ierr = PetscMalloc2(localSize,&dnnz,localSize,&onnz);CHKERRQ(ierr); 2360 2361 ierr = VecAssemblyEnd(vd_nz);CHKERRQ(ierr); 2362 ierr = VecAssemblyEnd(vo_nz);CHKERRQ(ierr); 2363 2364 ierr = VecGetArray(vd_nz,&vdnz);CHKERRQ(ierr); 2365 ierr = VecGetArray(vo_nz,&vonz);CHKERRQ(ierr); 2366 for (j=0; j<localSize; j++) { 2367 dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]); 2368 onnz[j] = (PetscInt)PetscRealPart(vonz[j]); 2369 } 2370 ierr = VecRestoreArray(vd_nz,&vdnz);CHKERRQ(ierr); 2371 ierr = VecRestoreArray(vo_nz,&vonz);CHKERRQ(ierr); 2372 ierr = VecDestroy(&vd_nz);CHKERRQ(ierr); 2373 ierr = VecDestroy(&vo_nz);CHKERRQ(ierr); 2374 2375 ierr = MatSeqAIJSetPreallocation(*J,0,dnnz);CHKERRQ(ierr); 2376 ierr = MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);CHKERRQ(ierr); 2377 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2378 2379 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 2380 2381 /* (2) Set matrix entries for edges */ 2382 /*----------------------------------*/ 2383 for (e=eStart; e<eEnd; e++) { 2384 /* Get row indices */ 2385 ierr = DMNetworkGetGlobalVecOffset(dm,e,ALL_COMPONENTS,&rstart);CHKERRQ(ierr); 2386 ierr = PetscSectionGetDof(network->DofSection,e,&nrows);CHKERRQ(ierr); 2387 if (nrows) { 2388 for (j=0; j<nrows; j++) rows[j] = j + rstart; 2389 2390 /* Set matrix entries for conntected vertices */ 2391 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 2392 for (v=0; v<2; v++) { 2393 ierr = DMNetworkGetGlobalVecOffset(dm,cone[v],ALL_COMPONENTS,&cstart);CHKERRQ(ierr); 2394 ierr = PetscSectionGetDof(network->DofSection,cone[v],&ncols);CHKERRQ(ierr); 2395 2396 if (network->Je) { 2397 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 2398 } else Juser = NULL; 2399 ierr = MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 2400 } 2401 2402 /* Set matrix entries for edge self */ 2403 cstart = rstart; 2404 if (network->Je) { 2405 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 2406 } else Juser = NULL; 2407 ierr = MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);CHKERRQ(ierr); 2408 } 2409 } 2410 2411 /* Set matrix entries for vertices */ 2412 /*---------------------------------*/ 2413 for (v=vStart; v<vEnd; v++) { 2414 /* Get row indices */ 2415 ierr = DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&rstart);CHKERRQ(ierr); 2416 ierr = PetscSectionGetDof(network->DofSection,v,&nrows);CHKERRQ(ierr); 2417 if (!nrows) continue; 2418 2419 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 2420 if (ghost) { 2421 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 2422 } else { 2423 rows_v = rows; 2424 } 2425 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 2426 2427 /* Get supporting edges and connected vertices */ 2428 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 2429 2430 for (e=0; e<nedges; e++) { 2431 /* Supporting edges */ 2432 ierr = DMNetworkGetGlobalVecOffset(dm,edges[e],ALL_COMPONENTS,&cstart);CHKERRQ(ierr); 2433 ierr = PetscSectionGetDof(network->DofSection,edges[e],&ncols);CHKERRQ(ierr); 2434 2435 if (network->Jv) { 2436 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 2437 } else Juser = NULL; 2438 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 2439 2440 /* Connected vertices */ 2441 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 2442 vc = (v == cone[0]) ? cone[1]:cone[0]; 2443 2444 ierr = DMNetworkGetGlobalVecOffset(dm,vc,ALL_COMPONENTS,&cstart);CHKERRQ(ierr); 2445 ierr = PetscSectionGetDof(network->DofSection,vc,&ncols);CHKERRQ(ierr); 2446 2447 if (network->Jv) { 2448 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 2449 } else Juser = NULL; 2450 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 2451 } 2452 2453 /* Set matrix entries for vertex self */ 2454 if (!ghost) { 2455 ierr = DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&cstart);CHKERRQ(ierr); 2456 if (network->Jv) { 2457 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 2458 } else Juser = NULL; 2459 ierr = MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);CHKERRQ(ierr); 2460 } 2461 if (ghost) { 2462 ierr = PetscFree(rows_v);CHKERRQ(ierr); 2463 } 2464 } 2465 ierr = PetscFree(rows);CHKERRQ(ierr); 2466 2467 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2468 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2469 2470 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 2471 PetscFunctionReturn(0); 2472 } 2473 2474 PetscErrorCode DMDestroy_Network(DM dm) 2475 { 2476 PetscErrorCode ierr; 2477 DM_Network *network = (DM_Network*)dm->data; 2478 PetscInt j; 2479 2480 PetscFunctionBegin; 2481 if (--network->refct > 0) PetscFunctionReturn(0); 2482 if (network->Je) { 2483 ierr = PetscFree(network->Je);CHKERRQ(ierr); 2484 } 2485 if (network->Jv) { 2486 ierr = PetscFree(network->Jvptr);CHKERRQ(ierr); 2487 ierr = PetscFree(network->Jv);CHKERRQ(ierr); 2488 } 2489 2490 ierr = ISLocalToGlobalMappingDestroy(&network->vertex.mapping);CHKERRQ(ierr); 2491 ierr = PetscSectionDestroy(&network->vertex.DofSection);CHKERRQ(ierr); 2492 ierr = PetscSectionDestroy(&network->vertex.GlobalDofSection);CHKERRQ(ierr); 2493 if (network->vltog) { 2494 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2495 } 2496 if (network->vertex.sf) { 2497 ierr = PetscSFDestroy(&network->vertex.sf);CHKERRQ(ierr); 2498 } 2499 /* edge */ 2500 ierr = ISLocalToGlobalMappingDestroy(&network->edge.mapping);CHKERRQ(ierr); 2501 ierr = PetscSectionDestroy(&network->edge.DofSection);CHKERRQ(ierr); 2502 ierr = PetscSectionDestroy(&network->edge.GlobalDofSection);CHKERRQ(ierr); 2503 if (network->edge.sf) { 2504 ierr = PetscSFDestroy(&network->edge.sf);CHKERRQ(ierr); 2505 } 2506 ierr = DMDestroy(&network->plex);CHKERRQ(ierr); 2507 ierr = PetscSectionDestroy(&network->DataSection);CHKERRQ(ierr); 2508 ierr = PetscSectionDestroy(&network->DofSection);CHKERRQ(ierr); 2509 2510 for (j=0; j<network->Nsvtx; j++) { 2511 ierr = PetscFree(network->svtx[j].sv);CHKERRQ(ierr); 2512 } 2513 if (network->svtx) {ierr = PetscFree(network->svtx);CHKERRQ(ierr);} 2514 2515 for (j=0; j<network->Nsubnet; j++) { 2516 ierr = PetscFree(network->subnet[j].edges);CHKERRQ(ierr); 2517 } 2518 if (network->subnetvtx) {ierr = PetscFree(network->subnetvtx);CHKERRQ(ierr);} 2519 2520 ierr = PetscTableDestroy(&network->svtable);CHKERRQ(ierr); 2521 ierr = PetscFree(network->subnet);CHKERRQ(ierr); 2522 ierr = PetscFree(network->componentdataarray);CHKERRQ(ierr); 2523 ierr = PetscFree2(network->header,network->cvalue);CHKERRQ(ierr); 2524 ierr = PetscFree(network);CHKERRQ(ierr); 2525 PetscFunctionReturn(0); 2526 } 2527 2528 PetscErrorCode DMView_Network(DM dm,PetscViewer viewer) 2529 { 2530 PetscErrorCode ierr; 2531 PetscBool iascii; 2532 PetscMPIInt rank; 2533 2534 PetscFunctionBegin; 2535 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 2536 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 2537 PetscValidHeaderSpecific(dm,DM_CLASSID, 1); 2538 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 2539 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2540 if (iascii) { 2541 const PetscInt *cone,*vtx,*edges; 2542 PetscInt vfrom,vto,i,j,nv,ne,ncv,p,nsubnet; 2543 DM_Network *network = (DM_Network*)dm->data; 2544 2545 nsubnet = network->Nsubnet; /* num of subnetworks */ 2546 if (!rank) { 2547 ierr = PetscPrintf(PETSC_COMM_SELF," NSubnets: %D; NEdges: %D; NVertices: %D; NSharedVertices: %D.\n",nsubnet,network->NEdges,network->NVertices,network->Nsvtx);CHKERRQ(ierr); 2548 } 2549 2550 ierr = DMNetworkGetSharedVertices(dm,&ncv,&vtx);CHKERRQ(ierr); 2551 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 2552 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nEdges: %D; nVertices: %D; nSharedVertices: %D\n",rank,network->nEdges,network->nVertices,ncv);CHKERRQ(ierr); 2553 2554 for (i=0; i<nsubnet; i++) { 2555 ierr = DMNetworkGetSubnetwork(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2556 if (ne) { 2557 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices(include shared vertices) %D\n",i,ne,nv);CHKERRQ(ierr); 2558 for (j=0; j<ne; j++) { 2559 p = edges[j]; 2560 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2561 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2562 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2563 ierr = DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);CHKERRQ(ierr); 2564 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D ----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2565 } 2566 } 2567 } 2568 2569 /* Shared vertices */ 2570 ierr = DMNetworkGetSharedVertices(dm,&ncv,&vtx);CHKERRQ(ierr); 2571 if (ncv) { 2572 SVtx *svtx = network->svtx; 2573 PetscInt gidx,svtx_idx,nvto,vfrom_net,vfrom_idx,*svto; 2574 PetscBool ghost; 2575 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " SharedVertices:\n");CHKERRQ(ierr); 2576 for (i=0; i<ncv; i++) { 2577 ierr = DMNetworkIsGhostVertex(dm,vtx[i],&ghost);CHKERRQ(ierr); 2578 if (ghost) continue; 2579 2580 ierr = DMNetworkGetGlobalVertexIndex(dm,vtx[i],&gidx);CHKERRQ(ierr); 2581 ierr = PetscTableFind(network->svtable,gidx+1,&svtx_idx);CHKERRQ(ierr); 2582 svtx_idx--; 2583 nvto = svtx[svtx_idx].n; 2584 2585 vfrom_net = svtx[svtx_idx].sv[0]; 2586 vfrom_idx = svtx[svtx_idx].sv[1]; 2587 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " svtx %D: global index %D, subnet[%D].%D ---->\n",i,gidx,vfrom_net,vfrom_idx);CHKERRQ(ierr); 2588 for (j=1; j<nvto; j++) { 2589 svto = svtx[svtx_idx].sv + 2*j; 2590 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " ----> subnet[%D].%D\n",svto[0],svto[1]);CHKERRQ(ierr); 2591 } 2592 } 2593 } 2594 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2595 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 2596 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name); 2597 PetscFunctionReturn(0); 2598 } 2599 2600 PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l) 2601 { 2602 PetscErrorCode ierr; 2603 DM_Network *network = (DM_Network*)dm->data; 2604 2605 PetscFunctionBegin; 2606 ierr = DMGlobalToLocalBegin(network->plex,g,mode,l);CHKERRQ(ierr); 2607 PetscFunctionReturn(0); 2608 } 2609 2610 PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l) 2611 { 2612 PetscErrorCode ierr; 2613 DM_Network *network = (DM_Network*)dm->data; 2614 2615 PetscFunctionBegin; 2616 ierr = DMGlobalToLocalEnd(network->plex,g,mode,l);CHKERRQ(ierr); 2617 PetscFunctionReturn(0); 2618 } 2619 2620 PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g) 2621 { 2622 PetscErrorCode ierr; 2623 DM_Network *network = (DM_Network*)dm->data; 2624 2625 PetscFunctionBegin; 2626 ierr = DMLocalToGlobalBegin(network->plex,l,mode,g);CHKERRQ(ierr); 2627 PetscFunctionReturn(0); 2628 } 2629 2630 PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g) 2631 { 2632 PetscErrorCode ierr; 2633 DM_Network *network = (DM_Network*)dm->data; 2634 2635 PetscFunctionBegin; 2636 ierr = DMLocalToGlobalEnd(network->plex,l,mode,g);CHKERRQ(ierr); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 /*@ 2641 DMNetworkGetVertexLocalToGlobalOrdering - Get vertex global index 2642 2643 Not collective 2644 2645 Input Parameters: 2646 + dm - the dm object 2647 - vloc - local vertex ordering, start from 0 2648 2649 Output Parameters: 2650 . vg - global vertex ordering, start from 0 2651 2652 Level: advanced 2653 2654 .seealso: DMNetworkSetVertexLocalToGlobalOrdering() 2655 @*/ 2656 PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg) 2657 { 2658 DM_Network *network = (DM_Network*)dm->data; 2659 PetscInt *vltog = network->vltog; 2660 2661 PetscFunctionBegin; 2662 if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first"); 2663 *vg = vltog[vloc]; 2664 PetscFunctionReturn(0); 2665 } 2666 2667 /*@ 2668 DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to global map 2669 2670 Collective 2671 2672 Input Parameters: 2673 . dm - the dm object 2674 2675 Level: advanced 2676 2677 .seealso: DMNetworkGetGlobalVertexIndex() 2678 @*/ 2679 PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm) 2680 { 2681 PetscErrorCode ierr; 2682 DM_Network *network=(DM_Network*)dm->data; 2683 MPI_Comm comm; 2684 PetscMPIInt rank,size,*displs=NULL,*recvcounts=NULL,remoterank; 2685 PetscBool ghost; 2686 PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx; 2687 const PetscSFNode *iremote; 2688 PetscSF vsf; 2689 Vec Vleaves,Vleaves_seq; 2690 VecScatter ctx; 2691 PetscScalar *varr,val; 2692 const PetscScalar *varr_read; 2693 2694 PetscFunctionBegin; 2695 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2696 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 2697 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 2698 2699 if (size == 1) { 2700 nroots = network->vEnd - network->vStart; 2701 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2702 for (i=0; i<nroots; i++) vltog[i] = i; 2703 network->vltog = vltog; 2704 PetscFunctionReturn(0); 2705 } 2706 2707 if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first"); 2708 if (network->vltog) { 2709 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2710 } 2711 2712 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 2713 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 2714 vsf = network->vertex.sf; 2715 2716 ierr = PetscMalloc3(size+1,&vrange,size,&displs,size,&recvcounts);CHKERRQ(ierr); 2717 ierr = PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);CHKERRQ(ierr); 2718 2719 for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;} 2720 2721 i = nroots - nleaves; /* local number of vertices, excluding ghosts */ 2722 vrange[0] = 0; 2723 ierr = MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);CHKERRMPI(ierr); 2724 for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];} 2725 2726 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2727 network->vltog = vltog; 2728 2729 /* Set vltog for non-ghost vertices */ 2730 k = 0; 2731 for (i=0; i<nroots; i++) { 2732 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2733 if (ghost) continue; 2734 vltog[i] = vrange[rank] + k++; 2735 } 2736 ierr = PetscFree3(vrange,displs,recvcounts);CHKERRQ(ierr); 2737 2738 /* Set vltog for ghost vertices */ 2739 /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */ 2740 ierr = VecCreate(comm,&Vleaves);CHKERRQ(ierr); 2741 ierr = VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);CHKERRQ(ierr); 2742 ierr = VecSetFromOptions(Vleaves);CHKERRQ(ierr); 2743 ierr = VecGetArray(Vleaves,&varr);CHKERRQ(ierr); 2744 for (i=0; i<nleaves; i++) { 2745 varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */ 2746 varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */ 2747 } 2748 ierr = VecRestoreArray(Vleaves,&varr);CHKERRQ(ierr); 2749 2750 /* (b) scatter local info to remote processes via VecScatter() */ 2751 ierr = VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);CHKERRQ(ierr); 2752 ierr = VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2753 ierr = VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2754 2755 /* (c) convert local indices to global indices in parallel vector Vleaves */ 2756 ierr = VecGetSize(Vleaves_seq,&N);CHKERRQ(ierr); 2757 ierr = VecGetArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2758 for (i=0; i<N; i+=2) { 2759 remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]); 2760 if (remoterank == rank) { 2761 k = i+1; /* row number */ 2762 lidx = (PetscInt)PetscRealPart(varr_read[i+1]); 2763 val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */ 2764 ierr = VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);CHKERRQ(ierr); 2765 } 2766 } 2767 ierr = VecRestoreArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2768 ierr = VecAssemblyBegin(Vleaves);CHKERRQ(ierr); 2769 ierr = VecAssemblyEnd(Vleaves);CHKERRQ(ierr); 2770 2771 /* (d) Set vltog for ghost vertices by copying local values of Vleaves */ 2772 ierr = VecGetArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2773 k = 0; 2774 for (i=0; i<nroots; i++) { 2775 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2776 if (!ghost) continue; 2777 vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++; 2778 } 2779 ierr = VecRestoreArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2780 2781 ierr = VecDestroy(&Vleaves);CHKERRQ(ierr); 2782 ierr = VecDestroy(&Vleaves_seq);CHKERRQ(ierr); 2783 ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr); 2784 PetscFunctionReturn(0); 2785 } 2786