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