1 #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/ 2 3 /*@ 4 DMNetworkGetPlex - Gets the Plex DM associated with this network DM 5 6 Not collective 7 8 Input Parameters: 9 + netdm - the dm object 10 - plexmdm - the plex dm object 11 12 Level: Advanced 13 14 .seealso: DMNetworkCreate() 15 @*/ 16 PetscErrorCode DMNetworkGetPlex(DM netdm, DM *plexdm) 17 { 18 DM_Network *network = (DM_Network*) netdm->data; 19 20 PetscFunctionBegin; 21 *plexdm = network->plex; 22 PetscFunctionReturn(0); 23 } 24 25 /*@ 26 DMNetworkGetSizes - Gets the the number of subnetworks and coupling subnetworks 27 28 Collective on dm 29 30 Input Parameters: 31 + dm - the dm object 32 . Nsubnet - global number of subnetworks 33 - NsubnetCouple - global number of coupling subnetworks 34 35 Level: Intermediate 36 37 .seealso: DMNetworkCreate() 38 @*/ 39 PetscErrorCode DMNetworkGetSizes(DM netdm, PetscInt *Nsubnet, PetscInt *Ncsubnet) 40 { 41 DM_Network *network = (DM_Network*) netdm->data; 42 43 PetscFunctionBegin; 44 *Nsubnet = network->nsubnet; 45 *Ncsubnet = network->ncsubnet; 46 PetscFunctionReturn(0); 47 } 48 49 /*@ 50 DMNetworkSetSizes - Sets the number of subnetworks, local and global vertices and edges for each subnetwork. 51 52 Collective on dm 53 54 Input Parameters: 55 + dm - the dm object 56 . Nsubnet - global number of subnetworks 57 . nV - number of local vertices for each subnetwork 58 . nE - number of local edges for each subnetwork 59 . NsubnetCouple - global number of coupling subnetworks 60 - nec - number of local edges for each coupling subnetwork 61 62 You cannot change the sizes once they have been set. nV, nE are arrays of length Nsubnet, and nec is array of length NsubnetCouple. 63 64 Level: intermediate 65 66 .seealso: DMNetworkCreate() 67 @*/ 68 PetscErrorCode DMNetworkSetSizes(DM dm,PetscInt Nsubnet,PetscInt nV[], PetscInt nE[],PetscInt NsubnetCouple,PetscInt nec[]) 69 { 70 PetscErrorCode ierr; 71 DM_Network *network = (DM_Network*) dm->data; 72 PetscInt a[2],b[2],i; 73 74 PetscFunctionBegin; 75 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 76 if (Nsubnet <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subnetworks %D cannot be less than 1",Nsubnet); 77 if (NsubnetCouple < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of coupling subnetworks %D cannot be less than 0",NsubnetCouple); 78 79 PetscValidLogicalCollectiveInt(dm,Nsubnet,2); 80 if (NsubnetCouple) PetscValidLogicalCollectiveInt(dm,NsubnetCouple,5); 81 if (network->nsubnet != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network"); 82 83 if (!nV || !nE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Local vertex size or edge size must be provided"); 84 85 network->nsubnet = Nsubnet + NsubnetCouple; 86 network->ncsubnet = NsubnetCouple; 87 ierr = PetscCalloc1(Nsubnet+NsubnetCouple,&network->subnet);CHKERRQ(ierr); 88 89 /* ---------------------------------------------------------- 90 p=v or e; P=V or E 91 subnet[0].pStart = 0 92 subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i]) 93 ----------------------------------------------------------------------- */ 94 for (i=0; i < Nsubnet; i++) { 95 /* Get global number of vertices and edges for subnet[i] */ 96 a[0] = nV[i]; a[1] = nE[i]; /* local number of vertices (excluding ghost) and edges */ 97 ierr = MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 98 network->subnet[i].Nvtx = b[0]; 99 network->subnet[i].Nedge = b[1]; 100 101 network->subnet[i].nvtx = nV[i]; /* local nvtx, without ghost */ 102 103 /* global subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */ 104 network->subnet[i].vStart = network->NVertices; 105 network->subnet[i].vEnd = network->subnet[i].vStart + network->subnet[i].Nvtx; 106 107 network->nVertices += nV[i]; 108 network->NVertices += network->subnet[i].Nvtx; 109 110 network->subnet[i].nedge = nE[i]; 111 network->subnet[i].eStart = network->nEdges; 112 network->subnet[i].eEnd = network->subnet[i].eStart + nE[i]; 113 network->nEdges += nE[i]; 114 network->NEdges += network->subnet[i].Nedge; 115 } 116 117 /* coupling subnetwork */ 118 for (; i < Nsubnet+NsubnetCouple; i++) { 119 /* Get global number of coupling edges for subnet[i] */ 120 ierr = MPIU_Allreduce(nec+(i-Nsubnet),&network->subnet[i].Nedge,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 121 122 network->subnet[i].nvtx = 0; /* We design coupling subnetwork such that it does not have its own vertices */ 123 network->subnet[i].vStart = network->nVertices; 124 network->subnet[i].vEnd = network->subnet[i].vStart; 125 126 network->subnet[i].nedge = nec[i-Nsubnet]; 127 network->subnet[i].eStart = network->nEdges; 128 network->subnet[i].eEnd = network->subnet[i].eStart + nec[i-Nsubnet]; 129 network->nEdges += nec[i-Nsubnet]; 130 network->NEdges += network->subnet[i].Nedge; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 /*@ 136 DMNetworkSetEdgeList - Sets the list of local edges (vertex connectivity) for the network 137 138 Logically collective on dm 139 140 Input Parameters: 141 + dm - the dm object 142 . edgelist - list of edges for each subnetwork 143 - edgelistCouple - list of edges for each coupling subnetwork 144 145 Notes: 146 There is no copy involved in this operation, only the pointer is referenced. The edgelist should 147 not be destroyed before the call to DMNetworkLayoutSetUp 148 149 Level: intermediate 150 151 Example usage: 152 Consider the following 2 separate networks and a coupling network: 153 154 .vb 155 network 0: v0 -> v1 -> v2 -> v3 156 network 1: v1 -> v2 -> v0 157 coupling network: network 1: v2 -> network 0: v0 158 .ve 159 160 The resulting input 161 edgelist[0] = [0 1 | 1 2 | 2 3]; 162 edgelist[1] = [1 2 | 2 0] 163 edgelistCouple[0] = [(network)1 (v)2 (network)0 (v)0]. 164 165 .seealso: DMNetworkCreate, DMNetworkSetSizes 166 @*/ 167 PetscErrorCode DMNetworkSetEdgeList(DM dm,PetscInt *edgelist[],PetscInt *edgelistCouple[]) 168 { 169 DM_Network *network = (DM_Network*) dm->data; 170 PetscInt i; 171 172 PetscFunctionBegin; 173 for (i=0; i < (network->nsubnet-network->ncsubnet); i++) network->subnet[i].edgelist = edgelist[i]; 174 if (network->ncsubnet) { 175 PetscInt j = 0; 176 if (!edgelistCouple) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Must provide edgelist_couple"); 177 while (i < network->nsubnet) network->subnet[i++].edgelist = edgelistCouple[j++]; 178 } 179 PetscFunctionReturn(0); 180 } 181 182 /*@ 183 DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network 184 185 Collective on dm 186 187 Input Parameters 188 . DM - the dmnetwork object 189 190 Notes: 191 This routine should be called after the network sizes and edgelists have been provided. It creates 192 the bare layout of the network and sets up the network to begin insertion of components. 193 194 All the components should be registered before calling this routine. 195 196 Level: intermediate 197 198 .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList 199 @*/ 200 PetscErrorCode DMNetworkLayoutSetUp(DM dm) 201 { 202 PetscErrorCode ierr; 203 DM_Network *network = (DM_Network*)dm->data; 204 PetscInt numCorners=2,spacedim=2,dim = 1; /* One dimensional network */ 205 PetscReal *vertexcoords=NULL; 206 PetscInt i,j,ctr,nsubnet,*eowners,np,*edges,*subnetvtx,vStart; 207 PetscInt k,netid,vid, *vidxlTog,*edgelist_couple=NULL; 208 const PetscInt *cone; 209 MPI_Comm comm; 210 PetscMPIInt size,rank; 211 212 PetscFunctionBegin; 213 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 214 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 215 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 216 217 /* Create the local edgelist for the network by concatenating local input edgelists of the subnetworks */ 218 ierr = PetscCalloc2(numCorners*network->nVertices,&vertexcoords,2*network->nEdges,&edges);CHKERRQ(ierr); 219 nsubnet = network->nsubnet - network->ncsubnet; 220 ctr = 0; 221 for (i=0; i < nsubnet; i++) { 222 for (j = 0; j < network->subnet[i].nedge; j++) { 223 edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j]; 224 edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1]; 225 ctr++; 226 } 227 } 228 229 /* Append local coupling edgelists of the subnetworks */ 230 i = nsubnet; /* netid of coupling subnet */ 231 nsubnet = network->nsubnet; 232 while (i < nsubnet) { 233 edgelist_couple = network->subnet[i].edgelist; 234 235 k = 0; 236 for (j = 0; j < network->subnet[i].nedge; j++) { 237 netid = edgelist_couple[k]; vid = edgelist_couple[k+1]; 238 edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2; 239 240 netid = edgelist_couple[k]; vid = edgelist_couple[k+1]; 241 edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2; 242 ctr++; 243 } 244 i++; 245 } 246 /* 247 if (rank == 0) { 248 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] edgelist:\n",rank); 249 for(i=0; i < network->nEdges; i++) { 250 ierr = PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edges[2*i],edges[2*i+1]);CHKERRQ(ierr); 251 printf("\n"); 252 } 253 } 254 */ 255 256 /* Create network->plex */ 257 #if defined(PETSC_USE_64BIT_INDICES) 258 { 259 int *edges64; 260 np = network->nEdges*numCorners; 261 ierr = PetscMalloc1(np,&edges64);CHKERRQ(ierr); 262 for (i=0; i<np; i++) edges64[i] = (int)edges[i]; 263 264 if (size == 1) { 265 ierr = DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const double*)vertexcoords,&network->plex);CHKERRQ(ierr); 266 } else { 267 ierr = DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);CHKERRQ(ierr); 268 } 269 ierr = PetscFree(edges64);CHKERRQ(ierr); 270 } 271 #else 272 if (size == 1) { 273 ierr = DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const double*)vertexcoords,&network->plex);CHKERRQ(ierr); 274 } else { 275 ierr = DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);CHKERRQ(ierr); 276 } 277 #endif 278 279 ierr = DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);CHKERRQ(ierr); 280 ierr = DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);CHKERRQ(ierr); 281 ierr = DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);CHKERRQ(ierr); 282 vStart = network->vStart; 283 284 ierr = PetscSectionCreate(comm,&network->DataSection);CHKERRQ(ierr); 285 ierr = PetscSectionCreate(comm,&network->DofSection);CHKERRQ(ierr); 286 ierr = PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);CHKERRQ(ierr); 287 ierr = PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);CHKERRQ(ierr); 288 289 network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType); 290 np = network->pEnd - network->pStart; 291 ierr = PetscCalloc2(np,&network->header,np,&network->cvalue);CHKERRQ(ierr); 292 293 /* Create vidxlTog: maps local vertex index to global index */ 294 np = network->vEnd - vStart; 295 ierr = PetscMalloc2(np,&vidxlTog,size+1,&eowners);CHKERRQ(ierr); 296 ctr = 0; 297 for (i=network->eStart; i<network->eEnd; i++) { 298 ierr = DMNetworkGetConnectedVertices(dm,i,&cone);CHKERRQ(ierr); 299 vidxlTog[cone[0] - vStart] = edges[2*ctr]; 300 vidxlTog[cone[1] - vStart] = edges[2*ctr+1]; 301 ctr++; 302 } 303 ierr = PetscFree2(vertexcoords,edges);CHKERRQ(ierr); 304 305 /* Create vertices and edges array for the subnetworks */ 306 for (j=0; j < network->nsubnet; j++) { 307 ierr = PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);CHKERRQ(ierr); 308 309 /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop. 310 These get updated when the vertices and edges are added. */ 311 network->subnet[j].nvtx = 0; 312 network->subnet[j].nedge = 0; 313 } 314 ierr = PetscCalloc1(np,&network->subnetvtx);CHKERRQ(ierr); 315 316 317 /* Get edge ownership */ 318 np = network->eEnd - network->eStart; 319 ierr = MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);CHKERRQ(ierr); 320 eowners[0] = 0; 321 for (i=2; i<=size; i++) eowners[i] += eowners[i-1]; 322 323 i = 0; j = 0; 324 while (i < np) { /* local edges, including coupling edges */ 325 network->header[i].index = i + eowners[rank]; /* Global edge index */ 326 327 if (j < network->nsubnet && i < network->subnet[j].eEnd) { 328 network->header[i].subnetid = j; /* Subnetwork id */ 329 network->subnet[j].edges[network->subnet[j].nedge++] = i; 330 331 network->header[i].ndata = 0; 332 ierr = PetscSectionAddDof(network->DataSection,i,network->dataheadersize);CHKERRQ(ierr); 333 network->header[i].offset[0] = 0; 334 i++; 335 } 336 if (i >= network->subnet[j].eEnd) j++; 337 } 338 339 /* Count network->subnet[*].nvtx */ 340 for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */ 341 k = vidxlTog[i-vStart]; 342 for (j=0; j < network->nsubnet; j++) { 343 if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) { 344 network->subnet[j].nvtx++; 345 break; 346 } 347 } 348 } 349 350 /* Set network->subnet[*].vertices on array network->subnetvtx */ 351 subnetvtx = network->subnetvtx; 352 for (j=0; j<network->nsubnet; j++) { 353 network->subnet[j].vertices = subnetvtx; 354 subnetvtx += network->subnet[j].nvtx; 355 network->subnet[j].nvtx = 0; 356 } 357 358 /* Set vertex array for the subnetworks */ 359 for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */ 360 network->header[i].index = vidxlTog[i-vStart]; /* Global vertex index */ 361 362 k = vidxlTog[i-vStart]; 363 for (j=0; j < network->nsubnet; j++) { 364 if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) { 365 network->header[i].subnetid = j; 366 network->subnet[j].vertices[network->subnet[j].nvtx++] = i; 367 break; 368 } 369 } 370 371 network->header[i].ndata = 0; 372 ierr = PetscSectionAddDof(network->DataSection,i,network->dataheadersize);CHKERRQ(ierr); 373 network->header[i].offset[0] = 0; 374 } 375 376 ierr = PetscFree2(vidxlTog,eowners);CHKERRQ(ierr); 377 PetscFunctionReturn(0); 378 } 379 380 /*@C 381 DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork 382 383 Input Parameters 384 + dm - the DM object 385 - id - the ID (integer) of the subnetwork 386 387 Output Parameters 388 + nv - number of vertices (local) 389 . ne - number of edges (local) 390 . vtx - local vertices for this subnetwork 391 - edge - local edges for this subnetwork 392 393 Notes: 394 Cannot call this routine before DMNetworkLayoutSetup() 395 396 Level: intermediate 397 398 .seealso: DMNetworkLayoutSetUp, DMNetworkCreate 399 @*/ 400 PetscErrorCode DMNetworkGetSubnetworkInfo(DM dm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge) 401 { 402 DM_Network *network = (DM_Network*)dm->data; 403 404 PetscFunctionBegin; 405 if (id >= network->nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of subnets %D",id,network->nsubnet); 406 *nv = network->subnet[id].nvtx; 407 *ne = network->subnet[id].nedge; 408 *vtx = network->subnet[id].vertices; 409 *edge = network->subnet[id].edges; 410 PetscFunctionReturn(0); 411 } 412 413 /*@C 414 DMNetworkGetSubnetworkCoupleInfo - Returns the info for the coupling subnetwork 415 416 Input Parameters 417 + dm - the DM object 418 - id - the ID (integer) of the coupling subnetwork 419 420 Output Parameters 421 + ne - number of edges (local) 422 - edge - local edges for this coupling subnetwork 423 424 Notes: 425 Cannot call this routine before DMNetworkLayoutSetup() 426 427 Level: intermediate 428 429 .seealso: DMNetworkGetSubnetworkInfo, DMNetworkLayoutSetUp, DMNetworkCreate 430 @*/ 431 PetscErrorCode DMNetworkGetSubnetworkCoupleInfo(DM dm,PetscInt id,PetscInt *ne,const PetscInt **edge) 432 { 433 DM_Network *net = (DM_Network*)dm->data; 434 PetscInt id1; 435 436 PetscFunctionBegin; 437 if (net->ncsubnet) { 438 if (id >= net->ncsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of coupling subnets %D",id,net->ncsubnet); 439 440 id1 = id + net->nsubnet - net->ncsubnet; 441 *ne = net->subnet[id1].nedge; 442 *edge = net->subnet[id1].edges; 443 } else { 444 *ne = 0; 445 *edge = NULL; 446 } 447 PetscFunctionReturn(0); 448 } 449 450 /*@C 451 DMNetworkRegisterComponent - Registers the network component 452 453 Logically collective on dm 454 455 Input Parameters 456 + dm - the network object 457 . name - the component name 458 - size - the storage size in bytes for this component data 459 460 Output Parameters 461 . key - an integer key that defines the component 462 463 Notes 464 This routine should be called by all processors before calling DMNetworkLayoutSetup(). 465 466 Level: intermediate 467 468 .seealso: DMNetworkLayoutSetUp, DMNetworkCreate 469 @*/ 470 PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key) 471 { 472 PetscErrorCode ierr; 473 DM_Network *network = (DM_Network*) dm->data; 474 DMNetworkComponent *component=&network->component[network->ncomponent]; 475 PetscBool flg=PETSC_FALSE; 476 PetscInt i; 477 478 PetscFunctionBegin; 479 for (i=0; i < network->ncomponent; i++) { 480 ierr = PetscStrcmp(component->name,name,&flg);CHKERRQ(ierr); 481 if (flg) { 482 *key = i; 483 PetscFunctionReturn(0); 484 } 485 } 486 if(network->ncomponent == MAX_COMPONENTS) { 487 SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS); 488 } 489 490 ierr = PetscStrcpy(component->name,name);CHKERRQ(ierr); 491 component->size = size/sizeof(DMNetworkComponentGenericDataType); 492 *key = network->ncomponent; 493 network->ncomponent++; 494 PetscFunctionReturn(0); 495 } 496 497 /*@ 498 DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices. 499 500 Not Collective 501 502 Input Parameters: 503 . dm - The DMNetwork object 504 505 Output Paramters: 506 + vStart - The first vertex point 507 - vEnd - One beyond the last vertex point 508 509 Level: intermediate 510 511 .seealso: DMNetworkGetEdgeRange 512 @*/ 513 PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd) 514 { 515 DM_Network *network = (DM_Network*)dm->data; 516 517 PetscFunctionBegin; 518 if (vStart) *vStart = network->vStart; 519 if (vEnd) *vEnd = network->vEnd; 520 PetscFunctionReturn(0); 521 } 522 523 /*@ 524 DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges. 525 526 Not Collective 527 528 Input Parameters: 529 . dm - The DMNetwork object 530 531 Output Paramters: 532 + eStart - The first edge point 533 - eEnd - One beyond the last edge point 534 535 Level: intermediate 536 537 .seealso: DMNetworkGetVertexRange 538 @*/ 539 PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd) 540 { 541 DM_Network *network = (DM_Network*)dm->data; 542 543 PetscFunctionBegin; 544 if (eStart) *eStart = network->eStart; 545 if (eEnd) *eEnd = network->eEnd; 546 PetscFunctionReturn(0); 547 } 548 549 /*@ 550 DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge. 551 552 Not Collective 553 554 Input Parameters: 555 + dm - DMNetwork object 556 - p - edge point 557 558 Output Paramters: 559 . index - user global numbering for the edge 560 561 Level: intermediate 562 563 .seealso: DMNetworkGetGlobalVertexIndex 564 @*/ 565 PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index) 566 { 567 PetscErrorCode ierr; 568 DM_Network *network = (DM_Network*)dm->data; 569 PetscInt offsetp; 570 DMNetworkComponentHeader header; 571 572 PetscFunctionBegin; 573 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 574 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 575 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp); 576 *index = header->index; 577 PetscFunctionReturn(0); 578 } 579 580 /*@ 581 DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex. 582 583 Not Collective 584 585 Input Parameters: 586 + dm - DMNetwork object 587 - p - vertex point 588 589 Output Paramters: 590 . index - user global numbering for the vertex 591 592 Level: intermediate 593 594 .seealso: DMNetworkGetGlobalEdgeIndex 595 @*/ 596 PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index) 597 { 598 PetscErrorCode ierr; 599 DM_Network *network = (DM_Network*)dm->data; 600 PetscInt offsetp; 601 DMNetworkComponentHeader header; 602 603 PetscFunctionBegin; 604 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 605 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 606 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp); 607 *index = header->index; 608 PetscFunctionReturn(0); 609 } 610 611 /* 612 DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the 613 component value from the component data array 614 615 Not Collective 616 617 Input Parameters: 618 + dm - The DMNetwork object 619 . p - vertex/edge point 620 - compnum - component number 621 622 Output Parameters: 623 + compkey - the key obtained when registering the component 624 - offset - offset into the component data array associated with the vertex/edge point 625 626 Notes: 627 Typical usage: 628 629 DMNetworkGetComponentDataArray(dm, &arr); 630 DMNetworkGetVertex/EdgeRange(dm,&Start,&End); 631 Loop over vertices or edges 632 DMNetworkGetNumComponents(dm,v,&numcomps); 633 Loop over numcomps 634 DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset); 635 compdata = (UserCompDataType)(arr+offset); 636 637 Level: intermediate 638 639 .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray, 640 */ 641 PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset) 642 { 643 PetscErrorCode ierr; 644 PetscInt offsetp; 645 DMNetworkComponentHeader header; 646 DM_Network *network = (DM_Network*)dm->data; 647 648 PetscFunctionBegin; 649 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 650 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp); 651 if (compkey) *compkey = header->key[compnum]; 652 if (offset) *offset = offsetp+network->dataheadersize+header->offset[compnum]; 653 PetscFunctionReturn(0); 654 } 655 656 /*@ 657 DMNetworkGetComponent - Returns the network component and its key 658 659 Not Collective 660 661 Input Parameters 662 + dm - DMNetwork object 663 . p - edge or vertex point 664 - compnum - component number 665 666 Output Parameters: 667 + compkey - the key set for this computing during registration 668 - component - the component data 669 670 Notes: 671 Typical usage: 672 673 DMNetworkGetVertex/EdgeRange(dm,&Start,&End); 674 Loop over vertices or edges 675 DMNetworkGetNumComponents(dm,v,&numcomps); 676 Loop over numcomps 677 DMNetworkGetComponent(dm,v,compnum,&key,&component); 678 679 Level: intermediate 680 681 .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset 682 @*/ 683 PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component) 684 { 685 PetscErrorCode ierr; 686 DM_Network *network = (DM_Network*)dm->data; 687 PetscInt offsetd = 0; 688 689 PetscFunctionBegin; 690 ierr = DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);CHKERRQ(ierr); 691 *component = network->componentdataarray+offsetd; 692 PetscFunctionReturn(0); 693 } 694 695 /*@ 696 DMNetworkAddComponent - Adds a network component at the given point (vertex/edge) 697 698 Not Collective 699 700 Input Parameters: 701 + dm - The DMNetwork object 702 . p - vertex/edge point 703 . componentkey - component key returned while registering the component 704 - compvalue - pointer to the data structure for the component 705 706 Level: intermediate 707 708 .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent 709 @*/ 710 PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue) 711 { 712 DM_Network *network = (DM_Network*)dm->data; 713 DMNetworkComponent *component = &network->component[componentkey]; 714 DMNetworkComponentHeader header = &network->header[p]; 715 DMNetworkComponentValue cvalue = &network->cvalue[p]; 716 PetscErrorCode ierr; 717 718 PetscFunctionBegin; 719 if (header->ndata == MAX_DATA_AT_POINT) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components at a point exceeds the max %D",MAX_DATA_AT_POINT); 720 721 header->size[header->ndata] = component->size; 722 ierr = PetscSectionAddDof(network->DataSection,p,component->size);CHKERRQ(ierr); 723 header->key[header->ndata] = componentkey; 724 if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1]; 725 726 cvalue->data[header->ndata] = (void*)compvalue; 727 header->ndata++; 728 PetscFunctionReturn(0); 729 } 730 731 /*@ 732 DMNetworkGetNumComponents - Get the number of components at a vertex/edge 733 734 Not Collective 735 736 Input Parameters: 737 + dm - The DMNetwork object 738 - p - vertex/edge point 739 740 Output Parameters: 741 . numcomponents - Number of components at the vertex/edge 742 743 Level: intermediate 744 745 .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent 746 @*/ 747 PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents) 748 { 749 PetscErrorCode ierr; 750 PetscInt offset; 751 DM_Network *network = (DM_Network*)dm->data; 752 753 PetscFunctionBegin; 754 ierr = PetscSectionGetOffset(network->DataSection,p,&offset);CHKERRQ(ierr); 755 *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata; 756 PetscFunctionReturn(0); 757 } 758 759 /*@ 760 DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector. 761 762 Not Collective 763 764 Input Parameters: 765 + dm - The DMNetwork object 766 - p - the edge/vertex point 767 768 Output Parameters: 769 . offset - the offset 770 771 Level: intermediate 772 773 .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector 774 @*/ 775 PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset) 776 { 777 PetscErrorCode ierr; 778 DM_Network *network = (DM_Network*)dm->data; 779 780 PetscFunctionBegin; 781 ierr = PetscSectionGetOffset(network->plex->localSection,p,offset);CHKERRQ(ierr); 782 PetscFunctionReturn(0); 783 } 784 785 /*@ 786 DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector. 787 788 Not Collective 789 790 Input Parameters: 791 + dm - The DMNetwork object 792 - p - the edge/vertex point 793 794 Output Parameters: 795 . offsetg - the offset 796 797 Level: intermediate 798 799 .seealso: DMNetworkGetVariableOffset, DMGetLocalVector 800 @*/ 801 PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg) 802 { 803 PetscErrorCode ierr; 804 DM_Network *network = (DM_Network*)dm->data; 805 806 PetscFunctionBegin; 807 ierr = PetscSectionGetOffset(network->plex->globalSection,p,offsetg);CHKERRQ(ierr); 808 if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */ 809 PetscFunctionReturn(0); 810 } 811 812 /*@ 813 DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector. 814 815 Not Collective 816 817 Input Parameters: 818 + dm - The DMNetwork object 819 - p - the edge point 820 821 Output Parameters: 822 . offset - the offset 823 824 Level: intermediate 825 826 .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector 827 @*/ 828 PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset) 829 { 830 PetscErrorCode ierr; 831 DM_Network *network = (DM_Network*)dm->data; 832 833 PetscFunctionBegin; 834 835 ierr = PetscSectionGetOffset(network->edge.DofSection,p,offset);CHKERRQ(ierr); 836 PetscFunctionReturn(0); 837 } 838 839 /*@ 840 DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector. 841 842 Not Collective 843 844 Input Parameters: 845 + dm - The DMNetwork object 846 - p - the vertex point 847 848 Output Parameters: 849 . offset - the offset 850 851 Level: intermediate 852 853 .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector 854 @*/ 855 PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset) 856 { 857 PetscErrorCode ierr; 858 DM_Network *network = (DM_Network*)dm->data; 859 860 PetscFunctionBegin; 861 862 p -= network->vStart; 863 864 ierr = PetscSectionGetOffset(network->vertex.DofSection,p,offset);CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867 /*@ 868 DMNetworkAddNumVariables - Add number of variables associated with a given point. 869 870 Not Collective 871 872 Input Parameters: 873 + dm - The DMNetworkObject 874 . p - the vertex/edge point 875 - nvar - number of additional variables 876 877 Level: intermediate 878 879 .seealso: DMNetworkSetNumVariables 880 @*/ 881 PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar) 882 { 883 PetscErrorCode ierr; 884 DM_Network *network = (DM_Network*)dm->data; 885 886 PetscFunctionBegin; 887 ierr = PetscSectionAddDof(network->DofSection,p,nvar);CHKERRQ(ierr); 888 PetscFunctionReturn(0); 889 } 890 891 /*@ 892 DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point. 893 894 Not Collective 895 896 Input Parameters: 897 + dm - The DMNetworkObject 898 - p - the vertex/edge point 899 900 Output Parameters: 901 . nvar - number of variables 902 903 Level: intermediate 904 905 .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables 906 @*/ 907 PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar) 908 { 909 PetscErrorCode ierr; 910 DM_Network *network = (DM_Network*)dm->data; 911 912 PetscFunctionBegin; 913 ierr = PetscSectionGetDof(network->DofSection,p,nvar);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 /*@ 918 DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point. 919 920 Not Collective 921 922 Input Parameters: 923 + dm - The DMNetworkObject 924 . p - the vertex/edge point 925 - nvar - number of variables 926 927 Level: intermediate 928 929 .seealso: DMNetworkAddNumVariables 930 @*/ 931 PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar) 932 { 933 PetscErrorCode ierr; 934 DM_Network *network = (DM_Network*)dm->data; 935 936 PetscFunctionBegin; 937 ierr = PetscSectionSetDof(network->DofSection,p,nvar);CHKERRQ(ierr); 938 PetscFunctionReturn(0); 939 } 940 941 /* Sets up the array that holds the data for all components and its associated section. This 942 function is called during DMSetUp() */ 943 PetscErrorCode DMNetworkComponentSetUp(DM dm) 944 { 945 PetscErrorCode ierr; 946 DM_Network *network = (DM_Network*)dm->data; 947 PetscInt arr_size,p,offset,offsetp,ncomp,i; 948 DMNetworkComponentHeader header; 949 DMNetworkComponentValue cvalue; 950 DMNetworkComponentGenericDataType *componentdataarray; 951 952 PetscFunctionBegin; 953 ierr = PetscSectionSetUp(network->DataSection);CHKERRQ(ierr); 954 ierr = PetscSectionGetStorageSize(network->DataSection,&arr_size);CHKERRQ(ierr); 955 ierr = PetscMalloc1(arr_size,&network->componentdataarray);CHKERRQ(ierr); 956 componentdataarray = network->componentdataarray; 957 for (p = network->pStart; p < network->pEnd; p++) { 958 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 959 /* Copy header */ 960 header = &network->header[p]; 961 ierr = PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));CHKERRQ(ierr); 962 /* Copy data */ 963 cvalue = &network->cvalue[p]; 964 ncomp = header->ndata; 965 for (i = 0; i < ncomp; i++) { 966 offset = offsetp + network->dataheadersize + header->offset[i]; 967 ierr = PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));CHKERRQ(ierr); 968 } 969 } 970 PetscFunctionReturn(0); 971 } 972 973 /* Sets up the section for dofs. This routine is called during DMSetUp() */ 974 PetscErrorCode DMNetworkVariablesSetUp(DM dm) 975 { 976 PetscErrorCode ierr; 977 DM_Network *network = (DM_Network*)dm->data; 978 979 PetscFunctionBegin; 980 ierr = PetscSectionSetUp(network->DofSection);CHKERRQ(ierr); 981 PetscFunctionReturn(0); 982 } 983 984 /*@C 985 DMNetworkGetComponentDataArray - Returns the component data array 986 987 Not Collective 988 989 Input Parameters: 990 . dm - The DMNetwork Object 991 992 Output Parameters: 993 . componentdataarray - array that holds data for all components 994 995 Level: intermediate 996 997 .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents 998 @*/ 999 PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray) 1000 { 1001 DM_Network *network = (DM_Network*)dm->data; 1002 1003 PetscFunctionBegin; 1004 *componentdataarray = network->componentdataarray; 1005 PetscFunctionReturn(0); 1006 } 1007 1008 /* Get a subsection from a range of points */ 1009 PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection) 1010 { 1011 PetscErrorCode ierr; 1012 PetscInt i, nvar; 1013 1014 PetscFunctionBegin; 1015 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);CHKERRQ(ierr); 1016 ierr = PetscSectionSetChart(*subsection, 0, pend - pstart);CHKERRQ(ierr); 1017 for (i = pstart; i < pend; i++) { 1018 ierr = PetscSectionGetDof(master,i,&nvar);CHKERRQ(ierr); 1019 ierr = PetscSectionSetDof(*subsection, i - pstart, nvar);CHKERRQ(ierr); 1020 } 1021 1022 ierr = PetscSectionSetUp(*subsection);CHKERRQ(ierr); 1023 PetscFunctionReturn(0); 1024 } 1025 1026 /* Create a submap of points with a GlobalToLocal structure */ 1027 PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map) 1028 { 1029 PetscErrorCode ierr; 1030 PetscInt i, *subpoints; 1031 1032 PetscFunctionBegin; 1033 /* Create index sets to map from "points" to "subpoints" */ 1034 ierr = PetscMalloc1(pend - pstart, &subpoints);CHKERRQ(ierr); 1035 for (i = pstart; i < pend; i++) { 1036 subpoints[i - pstart] = i; 1037 } 1038 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);CHKERRQ(ierr); 1039 ierr = PetscFree(subpoints);CHKERRQ(ierr); 1040 PetscFunctionReturn(0); 1041 } 1042 1043 /*@ 1044 DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute. 1045 1046 Collective 1047 1048 Input Parameters: 1049 . dm - The DMNetworkObject 1050 1051 Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are: 1052 1053 points = [0 1 2 3 4 5 6] 1054 1055 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]). 1056 1057 With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset. 1058 1059 Level: intermediate 1060 1061 @*/ 1062 PetscErrorCode DMNetworkAssembleGraphStructures(DM dm) 1063 { 1064 PetscErrorCode ierr; 1065 MPI_Comm comm; 1066 PetscMPIInt rank, size; 1067 DM_Network *network = (DM_Network*)dm->data; 1068 1069 PetscFunctionBegin; 1070 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1071 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1072 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1073 1074 /* Create maps for vertices and edges */ 1075 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 1076 ierr = DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);CHKERRQ(ierr); 1077 1078 /* Create local sub-sections */ 1079 ierr = DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);CHKERRQ(ierr); 1080 ierr = DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);CHKERRQ(ierr); 1081 1082 if (size > 1) { 1083 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 1084 1085 ierr = PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);CHKERRQ(ierr); 1086 ierr = PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);CHKERRQ(ierr); 1087 ierr = PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);CHKERRQ(ierr); 1088 } else { 1089 /* create structures for vertex */ 1090 ierr = PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);CHKERRQ(ierr); 1091 /* create structures for edge */ 1092 ierr = PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);CHKERRQ(ierr); 1093 } 1094 1095 /* Add viewers */ 1096 ierr = PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");CHKERRQ(ierr); 1097 ierr = PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");CHKERRQ(ierr); 1098 ierr = PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");CHKERRQ(ierr); 1099 ierr = PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");CHKERRQ(ierr); 1100 PetscFunctionReturn(0); 1101 } 1102 1103 /*@ 1104 DMNetworkDistribute - Distributes the network and moves associated component data. 1105 1106 Collective 1107 1108 Input Parameter: 1109 + DM - the DMNetwork object 1110 - overlap - The overlap of partitions, 0 is the default 1111 1112 Notes: 1113 Distributes the network with <overlap>-overlapping partitioning of the edges. 1114 1115 Level: intermediate 1116 1117 .seealso: DMNetworkCreate 1118 @*/ 1119 PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap) 1120 { 1121 MPI_Comm comm; 1122 PetscErrorCode ierr; 1123 PetscMPIInt size; 1124 DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data); 1125 DM_Network *newDMnetwork; 1126 PetscSF pointsf=NULL; 1127 DM newDM; 1128 PetscInt j,e,v,offset,*subnetvtx; 1129 PetscPartitioner part; 1130 DMNetworkComponentHeader header; 1131 1132 PetscFunctionBegin; 1133 ierr = PetscObjectGetComm((PetscObject)*dm,&comm);CHKERRQ(ierr); 1134 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1135 if (size == 1) PetscFunctionReturn(0); 1136 1137 ierr = DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);CHKERRQ(ierr); 1138 newDMnetwork = (DM_Network*)newDM->data; 1139 newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType); 1140 1141 /* Enable runtime options for petscpartitioner */ 1142 ierr = DMPlexGetPartitioner(oldDMnetwork->plex,&part);CHKERRQ(ierr); 1143 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 1144 1145 /* Distribute plex dm and dof section */ 1146 ierr = DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);CHKERRQ(ierr); 1147 1148 /* Distribute dof section */ 1149 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);CHKERRQ(ierr); 1150 ierr = PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);CHKERRQ(ierr); 1151 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);CHKERRQ(ierr); 1152 1153 /* Distribute data and associated section */ 1154 ierr = DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);CHKERRQ(ierr); 1155 1156 ierr = PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);CHKERRQ(ierr); 1157 ierr = DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);CHKERRQ(ierr); 1158 ierr = DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);CHKERRQ(ierr); 1159 newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart; 1160 newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart; 1161 newDMnetwork->NVertices = oldDMnetwork->NVertices; 1162 newDMnetwork->NEdges = oldDMnetwork->NEdges; 1163 1164 /* Set Dof section as the section for dm */ 1165 ierr = DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);CHKERRQ(ierr); 1166 ierr = DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);CHKERRQ(ierr); 1167 1168 /* Set up subnetwork info in the newDM */ 1169 newDMnetwork->nsubnet = oldDMnetwork->nsubnet; 1170 newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet; 1171 ierr = PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);CHKERRQ(ierr); 1172 /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already 1173 calculated in DMNetworkLayoutSetUp() 1174 */ 1175 for(j=0; j < newDMnetwork->nsubnet; j++) { 1176 newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx; 1177 newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge; 1178 } 1179 1180 for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) { 1181 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);CHKERRQ(ierr); 1182 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1183 newDMnetwork->subnet[header->subnetid].nedge++; 1184 } 1185 1186 for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) { 1187 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);CHKERRQ(ierr); 1188 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1189 newDMnetwork->subnet[header->subnetid].nvtx++; 1190 } 1191 1192 /* Now create the vertices and edge arrays for the subnetworks */ 1193 ierr = PetscCalloc1(newDMnetwork->vEnd-newDMnetwork->vStart,&newDMnetwork->subnetvtx);CHKERRQ(ierr); 1194 subnetvtx = newDMnetwork->subnetvtx; 1195 1196 for (j=0; j<newDMnetwork->nsubnet; j++) { 1197 ierr = PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);CHKERRQ(ierr); 1198 newDMnetwork->subnet[j].vertices = subnetvtx; 1199 subnetvtx += newDMnetwork->subnet[j].nvtx; 1200 1201 /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop. 1202 These get updated when the vertices and edges are added. */ 1203 newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0; 1204 } 1205 1206 /* Set the vertices and edges in each subnetwork */ 1207 for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) { 1208 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);CHKERRQ(ierr); 1209 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1210 newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e; 1211 } 1212 1213 for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) { 1214 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);CHKERRQ(ierr); 1215 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1216 newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v; 1217 } 1218 1219 newDM->setupcalled = (*dm)->setupcalled; 1220 newDMnetwork->distributecalled = PETSC_TRUE; 1221 1222 /* Destroy point SF */ 1223 ierr = PetscSFDestroy(&pointsf);CHKERRQ(ierr); 1224 1225 ierr = DMDestroy(dm);CHKERRQ(ierr); 1226 *dm = newDM; 1227 PetscFunctionReturn(0); 1228 } 1229 1230 /*@C 1231 PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering. 1232 1233 Input Parameters: 1234 + masterSF - the original SF structure 1235 - map - a ISLocalToGlobal mapping that contains the subset of points 1236 1237 Output Parameters: 1238 . subSF - a subset of the masterSF for the desired subset. 1239 */ 1240 1241 PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) { 1242 1243 PetscErrorCode ierr; 1244 PetscInt nroots, nleaves, *ilocal_sub; 1245 PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0; 1246 PetscInt *local_points, *remote_points; 1247 PetscSFNode *iremote_sub; 1248 const PetscInt *ilocal; 1249 const PetscSFNode *iremote; 1250 1251 PetscFunctionBegin; 1252 ierr = PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 1253 1254 /* Look for leaves that pertain to the subset of points. Get the local ordering */ 1255 ierr = PetscMalloc1(nleaves,&ilocal_map);CHKERRQ(ierr); 1256 ierr = ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);CHKERRQ(ierr); 1257 for (i = 0; i < nleaves; i++) { 1258 if (ilocal_map[i] != -1) nleaves_sub += 1; 1259 } 1260 /* Re-number ilocal with subset numbering. Need information from roots */ 1261 ierr = PetscMalloc2(nroots,&local_points,nroots,&remote_points);CHKERRQ(ierr); 1262 for (i = 0; i < nroots; i++) local_points[i] = i; 1263 ierr = ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);CHKERRQ(ierr); 1264 ierr = PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);CHKERRQ(ierr); 1265 ierr = PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);CHKERRQ(ierr); 1266 /* Fill up graph using local (that is, local to the subset) numbering. */ 1267 ierr = PetscMalloc1(nleaves_sub,&ilocal_sub);CHKERRQ(ierr); 1268 ierr = PetscMalloc1(nleaves_sub,&iremote_sub);CHKERRQ(ierr); 1269 nleaves_sub = 0; 1270 for (i = 0; i < nleaves; i++) { 1271 if (ilocal_map[i] != -1) { 1272 ilocal_sub[nleaves_sub] = ilocal_map[i]; 1273 iremote_sub[nleaves_sub].rank = iremote[i].rank; 1274 iremote_sub[nleaves_sub].index = remote_points[ilocal[i]]; 1275 nleaves_sub += 1; 1276 } 1277 } 1278 ierr = PetscFree2(local_points,remote_points);CHKERRQ(ierr); 1279 ierr = ISLocalToGlobalMappingGetSize(map,&nroots_sub);CHKERRQ(ierr); 1280 1281 /* Create new subSF */ 1282 ierr = PetscSFCreate(PETSC_COMM_WORLD,subSF);CHKERRQ(ierr); 1283 ierr = PetscSFSetFromOptions(*subSF);CHKERRQ(ierr); 1284 ierr = PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);CHKERRQ(ierr); 1285 ierr = PetscFree(ilocal_map);CHKERRQ(ierr); 1286 ierr = PetscFree(iremote_sub);CHKERRQ(ierr); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 /*@C 1291 DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point 1292 1293 Not Collective 1294 1295 Input Parameters: 1296 + dm - The DMNetwork object 1297 - p - the vertex point 1298 1299 Output Paramters: 1300 + nedges - number of edges connected to this vertex point 1301 - edges - List of edge points 1302 1303 Level: intermediate 1304 1305 Fortran Notes: 1306 Since it returns an array, this routine is only available in Fortran 90, and you must 1307 include petsc.h90 in your code. 1308 1309 .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices 1310 @*/ 1311 PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[]) 1312 { 1313 PetscErrorCode ierr; 1314 DM_Network *network = (DM_Network*)dm->data; 1315 1316 PetscFunctionBegin; 1317 ierr = DMPlexGetSupportSize(network->plex,vertex,nedges);CHKERRQ(ierr); 1318 ierr = DMPlexGetSupport(network->plex,vertex,edges);CHKERRQ(ierr); 1319 PetscFunctionReturn(0); 1320 } 1321 1322 /*@C 1323 DMNetworkGetConnectedVertices - Return the connected vertices for this edge point 1324 1325 Not Collective 1326 1327 Input Parameters: 1328 + dm - The DMNetwork object 1329 - p - the edge point 1330 1331 Output Paramters: 1332 . vertices - vertices connected to this edge 1333 1334 Level: intermediate 1335 1336 Fortran Notes: 1337 Since it returns an array, this routine is only available in Fortran 90, and you must 1338 include petsc.h90 in your code. 1339 1340 .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges 1341 @*/ 1342 PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[]) 1343 { 1344 PetscErrorCode ierr; 1345 DM_Network *network = (DM_Network*)dm->data; 1346 1347 PetscFunctionBegin; 1348 ierr = DMPlexGetCone(network->plex,edge,vertices);CHKERRQ(ierr); 1349 PetscFunctionReturn(0); 1350 } 1351 1352 /*@ 1353 DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex 1354 1355 Not Collective 1356 1357 Input Parameters: 1358 + dm - The DMNetwork object 1359 - p - the vertex point 1360 1361 Output Parameter: 1362 . isghost - TRUE if the vertex is a ghost point 1363 1364 Level: intermediate 1365 1366 .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange 1367 @*/ 1368 PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost) 1369 { 1370 PetscErrorCode ierr; 1371 DM_Network *network = (DM_Network*)dm->data; 1372 PetscInt offsetg; 1373 PetscSection sectiong; 1374 1375 PetscFunctionBegin; 1376 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 1377 *isghost = PETSC_FALSE; 1378 ierr = DMGetGlobalSection(network->plex,§iong);CHKERRQ(ierr); 1379 ierr = PetscSectionGetOffset(sectiong,p,&offsetg);CHKERRQ(ierr); 1380 if (offsetg < 0) *isghost = PETSC_TRUE; 1381 PetscFunctionReturn(0); 1382 } 1383 1384 PetscErrorCode DMSetUp_Network(DM dm) 1385 { 1386 PetscErrorCode ierr; 1387 DM_Network *network=(DM_Network*)dm->data; 1388 1389 PetscFunctionBegin; 1390 ierr = DMNetworkComponentSetUp(dm);CHKERRQ(ierr); 1391 ierr = DMNetworkVariablesSetUp(dm);CHKERRQ(ierr); 1392 1393 ierr = DMSetLocalSection(network->plex,network->DofSection);CHKERRQ(ierr); 1394 ierr = DMGetGlobalSection(network->plex,&network->GlobalDofSection);CHKERRQ(ierr); 1395 1396 dm->setupcalled = PETSC_TRUE; 1397 ierr = DMViewFromOptions(dm,NULL,"-dm_view");CHKERRQ(ierr); 1398 PetscFunctionReturn(0); 1399 } 1400 1401 /*@ 1402 DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices 1403 -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)? 1404 1405 Collective 1406 1407 Input Parameters: 1408 + dm - The DMNetwork object 1409 . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges 1410 - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices 1411 1412 Level: intermediate 1413 1414 @*/ 1415 PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg) 1416 { 1417 DM_Network *network=(DM_Network*)dm->data; 1418 PetscErrorCode ierr; 1419 PetscInt nVertices = network->nVertices; 1420 1421 PetscFunctionBegin; 1422 network->userEdgeJacobian = eflg; 1423 network->userVertexJacobian = vflg; 1424 1425 if (eflg && !network->Je) { 1426 ierr = PetscCalloc1(3*network->nEdges,&network->Je);CHKERRQ(ierr); 1427 } 1428 1429 if (vflg && !network->Jv && nVertices) { 1430 PetscInt i,*vptr,nedges,vStart=network->vStart; 1431 PetscInt nedges_total; 1432 const PetscInt *edges; 1433 1434 /* count nvertex_total */ 1435 nedges_total = 0; 1436 ierr = PetscMalloc1(nVertices+1,&vptr);CHKERRQ(ierr); 1437 1438 vptr[0] = 0; 1439 for (i=0; i<nVertices; i++) { 1440 ierr = DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);CHKERRQ(ierr); 1441 nedges_total += nedges; 1442 vptr[i+1] = vptr[i] + 2*nedges + 1; 1443 } 1444 1445 ierr = PetscCalloc1(2*nedges_total+nVertices,&network->Jv);CHKERRQ(ierr); 1446 network->Jvptr = vptr; 1447 } 1448 PetscFunctionReturn(0); 1449 } 1450 1451 /*@ 1452 DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network 1453 1454 Not Collective 1455 1456 Input Parameters: 1457 + dm - The DMNetwork object 1458 . p - the edge point 1459 - J - array (size = 3) of Jacobian submatrices for this edge point: 1460 J[0]: this edge 1461 J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices() 1462 1463 Level: intermediate 1464 1465 .seealso: DMNetworkVertexSetMatrix 1466 @*/ 1467 PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[]) 1468 { 1469 DM_Network *network=(DM_Network*)dm->data; 1470 1471 PetscFunctionBegin; 1472 if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix"); 1473 1474 if (J) { 1475 network->Je[3*p] = J[0]; 1476 network->Je[3*p+1] = J[1]; 1477 network->Je[3*p+2] = J[2]; 1478 } 1479 PetscFunctionReturn(0); 1480 } 1481 1482 /*@ 1483 DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network 1484 1485 Not Collective 1486 1487 Input Parameters: 1488 + dm - The DMNetwork object 1489 . p - the vertex point 1490 - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point: 1491 J[0]: this vertex 1492 J[1+2*i]: i-th supporting edge 1493 J[1+2*i+1]: i-th connected vertex 1494 1495 Level: intermediate 1496 1497 .seealso: DMNetworkEdgeSetMatrix 1498 @*/ 1499 PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[]) 1500 { 1501 PetscErrorCode ierr; 1502 DM_Network *network=(DM_Network*)dm->data; 1503 PetscInt i,*vptr,nedges,vStart=network->vStart; 1504 const PetscInt *edges; 1505 1506 PetscFunctionBegin; 1507 if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix"); 1508 1509 if (J) { 1510 vptr = network->Jvptr; 1511 network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */ 1512 1513 /* Set Jacobian for each supporting edge and connected vertex */ 1514 ierr = DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);CHKERRQ(ierr); 1515 for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i]; 1516 } 1517 PetscFunctionReturn(0); 1518 } 1519 1520 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 1521 { 1522 PetscErrorCode ierr; 1523 PetscInt j; 1524 PetscScalar val=(PetscScalar)ncols; 1525 1526 PetscFunctionBegin; 1527 if (!ghost) { 1528 for (j=0; j<nrows; j++) { 1529 ierr = VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1530 } 1531 } else { 1532 for (j=0; j<nrows; j++) { 1533 ierr = VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1534 } 1535 } 1536 PetscFunctionReturn(0); 1537 } 1538 1539 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 1540 { 1541 PetscErrorCode ierr; 1542 PetscInt j,ncols_u; 1543 PetscScalar val; 1544 1545 PetscFunctionBegin; 1546 if (!ghost) { 1547 for (j=0; j<nrows; j++) { 1548 ierr = MatGetRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1549 val = (PetscScalar)ncols_u; 1550 ierr = VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1551 ierr = MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1552 } 1553 } else { 1554 for (j=0; j<nrows; j++) { 1555 ierr = MatGetRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1556 val = (PetscScalar)ncols_u; 1557 ierr = VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1558 ierr = MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1559 } 1560 } 1561 PetscFunctionReturn(0); 1562 } 1563 1564 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 1565 { 1566 PetscErrorCode ierr; 1567 1568 PetscFunctionBegin; 1569 if (Ju) { 1570 ierr = MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);CHKERRQ(ierr); 1571 } else { 1572 ierr = MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);CHKERRQ(ierr); 1573 } 1574 PetscFunctionReturn(0); 1575 } 1576 1577 PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 1578 { 1579 PetscErrorCode ierr; 1580 PetscInt j,*cols; 1581 PetscScalar *zeros; 1582 1583 PetscFunctionBegin; 1584 ierr = PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);CHKERRQ(ierr); 1585 for (j=0; j<ncols; j++) cols[j] = j+ cstart; 1586 ierr = MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);CHKERRQ(ierr); 1587 ierr = PetscFree2(cols,zeros);CHKERRQ(ierr); 1588 PetscFunctionReturn(0); 1589 } 1590 1591 PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 1592 { 1593 PetscErrorCode ierr; 1594 PetscInt j,M,N,row,col,ncols_u; 1595 const PetscInt *cols; 1596 PetscScalar zero=0.0; 1597 1598 PetscFunctionBegin; 1599 ierr = MatGetSize(Ju,&M,&N);CHKERRQ(ierr); 1600 if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N); 1601 1602 for (row=0; row<nrows; row++) { 1603 ierr = MatGetRow(Ju,row,&ncols_u,&cols,NULL);CHKERRQ(ierr); 1604 for (j=0; j<ncols_u; j++) { 1605 col = cols[j] + cstart; 1606 ierr = MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);CHKERRQ(ierr); 1607 } 1608 ierr = MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);CHKERRQ(ierr); 1609 } 1610 PetscFunctionReturn(0); 1611 } 1612 1613 PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 1614 { 1615 PetscErrorCode ierr; 1616 1617 PetscFunctionBegin; 1618 if (Ju) { 1619 ierr = MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 1620 } else { 1621 ierr = MatSetDenseblock_private(nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 1622 } 1623 PetscFunctionReturn(0); 1624 } 1625 1626 /* 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. 1627 */ 1628 PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog) 1629 { 1630 PetscErrorCode ierr; 1631 PetscInt i,size,dof; 1632 PetscInt *glob2loc; 1633 1634 PetscFunctionBegin; 1635 ierr = PetscSectionGetStorageSize(localsec,&size);CHKERRQ(ierr); 1636 ierr = PetscMalloc1(size,&glob2loc);CHKERRQ(ierr); 1637 1638 for (i = 0; i < size; i++) { 1639 ierr = PetscSectionGetOffset(globalsec,i,&dof);CHKERRQ(ierr); 1640 dof = (dof >= 0) ? dof : -(dof + 1); 1641 glob2loc[i] = dof; 1642 } 1643 1644 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1645 #if 0 1646 ierr = PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1647 #endif 1648 PetscFunctionReturn(0); 1649 } 1650 1651 #include <petsc/private/matimpl.h> 1652 1653 PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J) 1654 { 1655 PetscErrorCode ierr; 1656 DM_Network *network = (DM_Network*)dm->data; 1657 PetscMPIInt rank, size; 1658 PetscInt eDof,vDof; 1659 Mat j11,j12,j21,j22,bA[2][2]; 1660 MPI_Comm comm; 1661 ISLocalToGlobalMapping eISMap,vISMap; 1662 1663 PetscFunctionBegin; 1664 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1665 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1666 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1667 1668 ierr = PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);CHKERRQ(ierr); 1669 ierr = PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);CHKERRQ(ierr); 1670 1671 ierr = MatCreate(comm, &j11);CHKERRQ(ierr); 1672 ierr = MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1673 ierr = MatSetType(j11, MATMPIAIJ);CHKERRQ(ierr); 1674 1675 ierr = MatCreate(comm, &j12);CHKERRQ(ierr); 1676 ierr = MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);CHKERRQ(ierr); 1677 ierr = MatSetType(j12, MATMPIAIJ);CHKERRQ(ierr); 1678 1679 ierr = MatCreate(comm, &j21);CHKERRQ(ierr); 1680 ierr = MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1681 ierr = MatSetType(j21, MATMPIAIJ);CHKERRQ(ierr); 1682 1683 ierr = MatCreate(comm, &j22);CHKERRQ(ierr); 1684 ierr = MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1685 ierr = MatSetType(j22, MATMPIAIJ);CHKERRQ(ierr); 1686 1687 bA[0][0] = j11; 1688 bA[0][1] = j12; 1689 bA[1][0] = j21; 1690 bA[1][1] = j22; 1691 1692 ierr = CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);CHKERRQ(ierr); 1693 ierr = CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);CHKERRQ(ierr); 1694 1695 ierr = MatSetLocalToGlobalMapping(j11,eISMap,eISMap);CHKERRQ(ierr); 1696 ierr = MatSetLocalToGlobalMapping(j12,eISMap,vISMap);CHKERRQ(ierr); 1697 ierr = MatSetLocalToGlobalMapping(j21,vISMap,eISMap);CHKERRQ(ierr); 1698 ierr = MatSetLocalToGlobalMapping(j22,vISMap,vISMap);CHKERRQ(ierr); 1699 1700 ierr = MatSetUp(j11);CHKERRQ(ierr); 1701 ierr = MatSetUp(j12);CHKERRQ(ierr); 1702 ierr = MatSetUp(j21);CHKERRQ(ierr); 1703 ierr = MatSetUp(j22);CHKERRQ(ierr); 1704 1705 ierr = MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);CHKERRQ(ierr); 1706 ierr = MatSetUp(*J);CHKERRQ(ierr); 1707 ierr = MatNestSetVecType(*J,VECNEST);CHKERRQ(ierr); 1708 ierr = MatDestroy(&j11);CHKERRQ(ierr); 1709 ierr = MatDestroy(&j12);CHKERRQ(ierr); 1710 ierr = MatDestroy(&j21);CHKERRQ(ierr); 1711 ierr = MatDestroy(&j22);CHKERRQ(ierr); 1712 1713 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1714 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1715 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1716 1717 /* Free structures */ 1718 ierr = ISLocalToGlobalMappingDestroy(&eISMap);CHKERRQ(ierr); 1719 ierr = ISLocalToGlobalMappingDestroy(&vISMap);CHKERRQ(ierr); 1720 PetscFunctionReturn(0); 1721 } 1722 1723 PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J) 1724 { 1725 PetscErrorCode ierr; 1726 DM_Network *network = (DM_Network*)dm->data; 1727 PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize; 1728 PetscInt cstart,ncols,j,e,v; 1729 PetscBool ghost,ghost_vc,ghost2,isNest; 1730 Mat Juser; 1731 PetscSection sectionGlobal; 1732 PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */ 1733 const PetscInt *edges,*cone; 1734 MPI_Comm comm; 1735 MatType mtype; 1736 Vec vd_nz,vo_nz; 1737 PetscInt *dnnz,*onnz; 1738 PetscScalar *vdnz,*vonz; 1739 1740 PetscFunctionBegin; 1741 mtype = dm->mattype; 1742 ierr = PetscStrcmp(mtype,MATNEST,&isNest);CHKERRQ(ierr); 1743 if (isNest) { 1744 ierr = DMCreateMatrix_Network_Nest(dm,J);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 if (!network->userEdgeJacobian && !network->userVertexJacobian) { 1749 /* user does not provide Jacobian blocks */ 1750 ierr = DMCreateMatrix_Plex(network->plex,J);CHKERRQ(ierr); 1751 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 1752 PetscFunctionReturn(0); 1753 } 1754 1755 ierr = MatCreate(PetscObjectComm((PetscObject)dm),J);CHKERRQ(ierr); 1756 ierr = DMGetGlobalSection(network->plex,§ionGlobal);CHKERRQ(ierr); 1757 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);CHKERRQ(ierr); 1758 ierr = MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1759 1760 ierr = MatSetType(*J,MATAIJ);CHKERRQ(ierr); 1761 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 1762 1763 /* (1) Set matrix preallocation */ 1764 /*------------------------------*/ 1765 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1766 ierr = VecCreate(comm,&vd_nz);CHKERRQ(ierr); 1767 ierr = VecSetSizes(vd_nz,localSize,PETSC_DECIDE);CHKERRQ(ierr); 1768 ierr = VecSetFromOptions(vd_nz);CHKERRQ(ierr); 1769 ierr = VecSet(vd_nz,0.0);CHKERRQ(ierr); 1770 ierr = VecDuplicate(vd_nz,&vo_nz);CHKERRQ(ierr); 1771 1772 /* Set preallocation for edges */ 1773 /*-----------------------------*/ 1774 ierr = DMNetworkGetEdgeRange(dm,&eStart,&eEnd);CHKERRQ(ierr); 1775 1776 ierr = PetscMalloc1(localSize,&rows);CHKERRQ(ierr); 1777 for (e=eStart; e<eEnd; e++) { 1778 /* Get row indices */ 1779 ierr = DMNetworkGetVariableGlobalOffset(dm,e,&rstart);CHKERRQ(ierr); 1780 ierr = DMNetworkGetNumVariables(dm,e,&nrows);CHKERRQ(ierr); 1781 if (nrows) { 1782 for (j=0; j<nrows; j++) rows[j] = j + rstart; 1783 1784 /* Set preallocation for conntected vertices */ 1785 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 1786 for (v=0; v<2; v++) { 1787 ierr = DMNetworkGetNumVariables(dm,cone[v],&ncols);CHKERRQ(ierr); 1788 1789 if (network->Je) { 1790 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 1791 } else Juser = NULL; 1792 ierr = DMNetworkIsGhostVertex(dm,cone[v],&ghost);CHKERRQ(ierr); 1793 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 1794 } 1795 1796 /* Set preallocation for edge self */ 1797 cstart = rstart; 1798 if (network->Je) { 1799 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 1800 } else Juser = NULL; 1801 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 1802 } 1803 } 1804 1805 /* Set preallocation for vertices */ 1806 /*--------------------------------*/ 1807 ierr = DMNetworkGetVertexRange(dm,&vStart,&vEnd);CHKERRQ(ierr); 1808 if (vEnd - vStart) vptr = network->Jvptr; 1809 1810 for (v=vStart; v<vEnd; v++) { 1811 /* Get row indices */ 1812 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&rstart);CHKERRQ(ierr); 1813 ierr = DMNetworkGetNumVariables(dm,v,&nrows);CHKERRQ(ierr); 1814 if (!nrows) continue; 1815 1816 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1817 if (ghost) { 1818 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 1819 } else { 1820 rows_v = rows; 1821 } 1822 1823 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 1824 1825 /* Get supporting edges and connected vertices */ 1826 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 1827 1828 for (e=0; e<nedges; e++) { 1829 /* Supporting edges */ 1830 ierr = DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);CHKERRQ(ierr); 1831 ierr = DMNetworkGetNumVariables(dm,edges[e],&ncols);CHKERRQ(ierr); 1832 1833 if (network->Jv) { 1834 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 1835 } else Juser = NULL; 1836 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 1837 1838 /* Connected vertices */ 1839 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 1840 vc = (v == cone[0]) ? cone[1]:cone[0]; 1841 ierr = DMNetworkIsGhostVertex(dm,vc,&ghost_vc);CHKERRQ(ierr); 1842 1843 ierr = DMNetworkGetNumVariables(dm,vc,&ncols);CHKERRQ(ierr); 1844 1845 if (network->Jv) { 1846 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 1847 } else Juser = NULL; 1848 if (ghost_vc||ghost) { 1849 ghost2 = PETSC_TRUE; 1850 } else { 1851 ghost2 = PETSC_FALSE; 1852 } 1853 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);CHKERRQ(ierr); 1854 } 1855 1856 /* Set preallocation for vertex self */ 1857 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1858 if (!ghost) { 1859 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&cstart);CHKERRQ(ierr); 1860 if (network->Jv) { 1861 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 1862 } else Juser = NULL; 1863 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 1864 } 1865 if (ghost) { 1866 ierr = PetscFree(rows_v);CHKERRQ(ierr); 1867 } 1868 } 1869 1870 ierr = VecAssemblyBegin(vd_nz);CHKERRQ(ierr); 1871 ierr = VecAssemblyBegin(vo_nz);CHKERRQ(ierr); 1872 1873 ierr = PetscMalloc2(localSize,&dnnz,localSize,&onnz);CHKERRQ(ierr); 1874 1875 ierr = VecAssemblyEnd(vd_nz);CHKERRQ(ierr); 1876 ierr = VecAssemblyEnd(vo_nz);CHKERRQ(ierr); 1877 1878 ierr = VecGetArray(vd_nz,&vdnz);CHKERRQ(ierr); 1879 ierr = VecGetArray(vo_nz,&vonz);CHKERRQ(ierr); 1880 for (j=0; j<localSize; j++) { 1881 dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]); 1882 onnz[j] = (PetscInt)PetscRealPart(vonz[j]); 1883 } 1884 ierr = VecRestoreArray(vd_nz,&vdnz);CHKERRQ(ierr); 1885 ierr = VecRestoreArray(vo_nz,&vonz);CHKERRQ(ierr); 1886 ierr = VecDestroy(&vd_nz);CHKERRQ(ierr); 1887 ierr = VecDestroy(&vo_nz);CHKERRQ(ierr); 1888 1889 ierr = MatSeqAIJSetPreallocation(*J,0,dnnz);CHKERRQ(ierr); 1890 ierr = MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);CHKERRQ(ierr); 1891 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1892 1893 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 1894 1895 /* (2) Set matrix entries for edges */ 1896 /*----------------------------------*/ 1897 for (e=eStart; e<eEnd; e++) { 1898 /* Get row indices */ 1899 ierr = DMNetworkGetVariableGlobalOffset(dm,e,&rstart);CHKERRQ(ierr); 1900 ierr = DMNetworkGetNumVariables(dm,e,&nrows);CHKERRQ(ierr); 1901 if (nrows) { 1902 for (j=0; j<nrows; j++) rows[j] = j + rstart; 1903 1904 /* Set matrix entries for conntected vertices */ 1905 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 1906 for (v=0; v<2; v++) { 1907 ierr = DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);CHKERRQ(ierr); 1908 ierr = DMNetworkGetNumVariables(dm,cone[v],&ncols);CHKERRQ(ierr); 1909 1910 if (network->Je) { 1911 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 1912 } else Juser = NULL; 1913 ierr = MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 1914 } 1915 1916 /* Set matrix entries for edge self */ 1917 cstart = rstart; 1918 if (network->Je) { 1919 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 1920 } else Juser = NULL; 1921 ierr = MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);CHKERRQ(ierr); 1922 } 1923 } 1924 1925 /* Set matrix entries for vertices */ 1926 /*---------------------------------*/ 1927 for (v=vStart; v<vEnd; v++) { 1928 /* Get row indices */ 1929 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&rstart);CHKERRQ(ierr); 1930 ierr = DMNetworkGetNumVariables(dm,v,&nrows);CHKERRQ(ierr); 1931 if (!nrows) continue; 1932 1933 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1934 if (ghost) { 1935 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 1936 } else { 1937 rows_v = rows; 1938 } 1939 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 1940 1941 /* Get supporting edges and connected vertices */ 1942 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 1943 1944 for (e=0; e<nedges; e++) { 1945 /* Supporting edges */ 1946 ierr = DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);CHKERRQ(ierr); 1947 ierr = DMNetworkGetNumVariables(dm,edges[e],&ncols);CHKERRQ(ierr); 1948 1949 if (network->Jv) { 1950 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 1951 } else Juser = NULL; 1952 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 1953 1954 /* Connected vertices */ 1955 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 1956 vc = (v == cone[0]) ? cone[1]:cone[0]; 1957 1958 ierr = DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);CHKERRQ(ierr); 1959 ierr = DMNetworkGetNumVariables(dm,vc,&ncols);CHKERRQ(ierr); 1960 1961 if (network->Jv) { 1962 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 1963 } else Juser = NULL; 1964 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 1965 } 1966 1967 /* Set matrix entries for vertex self */ 1968 if (!ghost) { 1969 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&cstart);CHKERRQ(ierr); 1970 if (network->Jv) { 1971 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 1972 } else Juser = NULL; 1973 ierr = MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);CHKERRQ(ierr); 1974 } 1975 if (ghost) { 1976 ierr = PetscFree(rows_v);CHKERRQ(ierr); 1977 } 1978 } 1979 ierr = PetscFree(rows);CHKERRQ(ierr); 1980 1981 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1982 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1983 1984 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 1985 PetscFunctionReturn(0); 1986 } 1987 1988 PetscErrorCode DMDestroy_Network(DM dm) 1989 { 1990 PetscErrorCode ierr; 1991 DM_Network *network = (DM_Network*)dm->data; 1992 PetscInt j; 1993 1994 PetscFunctionBegin; 1995 if (--network->refct > 0) PetscFunctionReturn(0); 1996 if (network->Je) { 1997 ierr = PetscFree(network->Je);CHKERRQ(ierr); 1998 } 1999 if (network->Jv) { 2000 ierr = PetscFree(network->Jvptr);CHKERRQ(ierr); 2001 ierr = PetscFree(network->Jv);CHKERRQ(ierr); 2002 } 2003 2004 ierr = ISLocalToGlobalMappingDestroy(&network->vertex.mapping);CHKERRQ(ierr); 2005 ierr = PetscSectionDestroy(&network->vertex.DofSection);CHKERRQ(ierr); 2006 ierr = PetscSectionDestroy(&network->vertex.GlobalDofSection);CHKERRQ(ierr); 2007 if (network->vltog) { 2008 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2009 } 2010 if (network->vertex.sf) { 2011 ierr = PetscSFDestroy(&network->vertex.sf);CHKERRQ(ierr); 2012 } 2013 /* edge */ 2014 ierr = ISLocalToGlobalMappingDestroy(&network->edge.mapping);CHKERRQ(ierr); 2015 ierr = PetscSectionDestroy(&network->edge.DofSection);CHKERRQ(ierr); 2016 ierr = PetscSectionDestroy(&network->edge.GlobalDofSection);CHKERRQ(ierr); 2017 if (network->edge.sf) { 2018 ierr = PetscSFDestroy(&network->edge.sf);CHKERRQ(ierr); 2019 } 2020 ierr = DMDestroy(&network->plex);CHKERRQ(ierr); 2021 ierr = PetscSectionDestroy(&network->DataSection);CHKERRQ(ierr); 2022 ierr = PetscSectionDestroy(&network->DofSection);CHKERRQ(ierr); 2023 2024 for(j=0; j<network->nsubnet; j++) { 2025 ierr = PetscFree(network->subnet[j].edges);CHKERRQ(ierr); 2026 } 2027 ierr = PetscFree(network->subnetvtx);CHKERRQ(ierr); 2028 2029 ierr = PetscFree(network->subnet);CHKERRQ(ierr); 2030 ierr = PetscFree(network->componentdataarray);CHKERRQ(ierr); 2031 ierr = PetscFree2(network->header,network->cvalue);CHKERRQ(ierr); 2032 ierr = PetscFree(network);CHKERRQ(ierr); 2033 PetscFunctionReturn(0); 2034 } 2035 2036 PetscErrorCode DMView_Network(DM dm,PetscViewer viewer) 2037 { 2038 PetscErrorCode ierr; 2039 DM_Network *network = (DM_Network*)dm->data; 2040 PetscBool iascii; 2041 PetscMPIInt rank; 2042 PetscInt p,nsubnet; 2043 2044 PetscFunctionBegin; 2045 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 2046 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 2047 PetscValidHeaderSpecific(dm,DM_CLASSID, 1); 2048 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 2049 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2050 if (iascii) { 2051 const PetscInt *cone,*vtx,*edges; 2052 PetscInt vfrom,vto,i,j,nv,ne; 2053 2054 nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */ 2055 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 2056 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);CHKERRQ(ierr); 2057 2058 for (i=0; i<nsubnet; i++) { 2059 ierr = DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2060 if (ne) { 2061 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);CHKERRQ(ierr); 2062 for (j=0; j<ne; j++) { 2063 p = edges[j]; 2064 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2065 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2066 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2067 ierr = DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);CHKERRQ(ierr); 2068 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2069 } 2070 } 2071 } 2072 /* Coupling subnets */ 2073 nsubnet = network->nsubnet; 2074 for (; i<nsubnet; i++) { 2075 ierr = DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2076 if (ne) { 2077 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);CHKERRQ(ierr); 2078 for (j=0; j<ne; j++) { 2079 p = edges[j]; 2080 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2081 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2082 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2083 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2084 } 2085 } 2086 } 2087 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2088 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 2089 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l) 2094 { 2095 PetscErrorCode ierr; 2096 DM_Network *network = (DM_Network*)dm->data; 2097 2098 PetscFunctionBegin; 2099 ierr = DMGlobalToLocalBegin(network->plex,g,mode,l);CHKERRQ(ierr); 2100 PetscFunctionReturn(0); 2101 } 2102 2103 PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l) 2104 { 2105 PetscErrorCode ierr; 2106 DM_Network *network = (DM_Network*)dm->data; 2107 2108 PetscFunctionBegin; 2109 ierr = DMGlobalToLocalEnd(network->plex,g,mode,l);CHKERRQ(ierr); 2110 PetscFunctionReturn(0); 2111 } 2112 2113 PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g) 2114 { 2115 PetscErrorCode ierr; 2116 DM_Network *network = (DM_Network*)dm->data; 2117 2118 PetscFunctionBegin; 2119 ierr = DMLocalToGlobalBegin(network->plex,l,mode,g);CHKERRQ(ierr); 2120 PetscFunctionReturn(0); 2121 } 2122 2123 PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g) 2124 { 2125 PetscErrorCode ierr; 2126 DM_Network *network = (DM_Network*)dm->data; 2127 2128 PetscFunctionBegin; 2129 ierr = DMLocalToGlobalEnd(network->plex,l,mode,g);CHKERRQ(ierr); 2130 PetscFunctionReturn(0); 2131 } 2132 2133 /*@ 2134 DMNetworkGetVertexLocalToGlobalOrdering - Get vertex globle index 2135 2136 Not collective 2137 2138 Input Parameters 2139 + dm - the dm object 2140 - vloc - local vertex ordering, start from 0 2141 2142 Output Parameters 2143 + vg - global vertex ordering, start from 0 2144 2145 Level: Advanced 2146 2147 .seealso: DMNetworkSetVertexLocalToGlobalOrdering() 2148 @*/ 2149 PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg) 2150 { 2151 DM_Network *network = (DM_Network*)dm->data; 2152 PetscInt *vltog = network->vltog; 2153 2154 PetscFunctionBegin; 2155 if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first"); 2156 *vg = vltog[vloc]; 2157 PetscFunctionReturn(0); 2158 } 2159 2160 /*@ 2161 DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to globle map 2162 2163 Collective 2164 2165 Input Parameters: 2166 + dm - the dm object 2167 2168 Level: Advanced 2169 2170 .seealso: DMNetworkGetGlobalVertexIndex() 2171 @*/ 2172 PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm) 2173 { 2174 PetscErrorCode ierr; 2175 DM_Network *network=(DM_Network*)dm->data; 2176 MPI_Comm comm; 2177 PetscMPIInt rank,size,*displs,*recvcounts,remoterank; 2178 PetscBool ghost; 2179 PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx; 2180 const PetscSFNode *iremote; 2181 PetscSF vsf; 2182 Vec Vleaves,Vleaves_seq; 2183 VecScatter ctx; 2184 PetscScalar *varr,val; 2185 const PetscScalar *varr_read; 2186 2187 PetscFunctionBegin; 2188 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2189 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2190 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 2191 2192 if (size == 1) { 2193 nroots = network->vEnd - network->vStart; 2194 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2195 for (i=0; i<nroots; i++) vltog[i] = i; 2196 network->vltog = vltog; 2197 PetscFunctionReturn(0); 2198 } 2199 2200 if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first"); 2201 if (network->vltog) { 2202 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2203 } 2204 2205 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 2206 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 2207 vsf = network->vertex.sf; 2208 2209 ierr = PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);CHKERRQ(ierr); 2210 ierr = PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);CHKERRQ(ierr); 2211 2212 for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;} 2213 2214 i = nroots - nleaves; /* local number of vertices, excluding ghosts */ 2215 vrange[0] = 0; 2216 ierr = MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);CHKERRQ(ierr); 2217 for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];} 2218 2219 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2220 network->vltog = vltog; 2221 2222 /* Set vltog for non-ghost vertices */ 2223 k = 0; 2224 for (i=0; i<nroots; i++) { 2225 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2226 if (ghost) continue; 2227 vltog[i] = vrange[rank] + k++; 2228 } 2229 ierr = PetscFree3(vrange,displs,recvcounts);CHKERRQ(ierr); 2230 2231 /* Set vltog for ghost vertices */ 2232 /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */ 2233 ierr = VecCreate(comm,&Vleaves);CHKERRQ(ierr); 2234 ierr = VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);CHKERRQ(ierr); 2235 ierr = VecSetFromOptions(Vleaves);CHKERRQ(ierr); 2236 ierr = VecGetArray(Vleaves,&varr);CHKERRQ(ierr); 2237 for (i=0; i<nleaves; i++) { 2238 varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */ 2239 varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */ 2240 } 2241 ierr = VecRestoreArray(Vleaves,&varr);CHKERRQ(ierr); 2242 2243 /* (b) scatter local info to remote processes via VecScatter() */ 2244 ierr = VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);CHKERRQ(ierr); 2245 ierr = VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2246 ierr = VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2247 2248 /* (c) convert local indices to global indices in parallel vector Vleaves */ 2249 ierr = VecGetSize(Vleaves_seq,&N);CHKERRQ(ierr); 2250 ierr = VecGetArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2251 for (i=0; i<N; i+=2) { 2252 remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]); 2253 if (remoterank == rank) { 2254 k = i+1; /* row number */ 2255 lidx = (PetscInt)PetscRealPart(varr_read[i+1]); 2256 val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */ 2257 ierr = VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);CHKERRQ(ierr); 2258 } 2259 } 2260 ierr = VecRestoreArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2261 ierr = VecAssemblyBegin(Vleaves);CHKERRQ(ierr); 2262 ierr = VecAssemblyEnd(Vleaves);CHKERRQ(ierr); 2263 2264 /* (d) Set vltog for ghost vertices by copying local values of Vleaves */ 2265 ierr = VecGetArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2266 k = 0; 2267 for (i=0; i<nroots; i++) { 2268 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2269 if (!ghost) continue; 2270 vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++; 2271 } 2272 ierr = VecRestoreArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2273 2274 ierr = VecDestroy(&Vleaves);CHKERRQ(ierr); 2275 ierr = VecDestroy(&Vleaves_seq);CHKERRQ(ierr); 2276 ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr); 2277 PetscFunctionReturn(0); 2278 } 2279