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