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