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