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 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: beginner 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 Parameters: 508 + vStart - The first vertex point 509 - vEnd - One beyond the last vertex point 510 511 Level: beginner 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 Parameters: 534 + eStart - The first edge point 535 - eEnd - One beyond the last edge point 536 537 Level: beginner 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 Parameters: 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 Parameters: 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: beginner 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: beginner 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: beginner 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 ierr = DMNetworkAddNumVariables(dm,p,nvar);CHKERRQ(ierr); 757 header->nvar[compnum] = nvar; 758 if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1]; 759 PetscFunctionReturn(0); 760 } 761 762 /*@ 763 DMNetworkGetNumComponents - Get the number of components at a vertex/edge 764 765 Not Collective 766 767 Input Parameters: 768 + dm - The DMNetwork object 769 - p - vertex/edge point 770 771 Output Parameters: 772 . numcomponents - Number of components at the vertex/edge 773 774 Level: beginner 775 776 .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent 777 @*/ 778 PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents) 779 { 780 PetscErrorCode ierr; 781 PetscInt offset; 782 DM_Network *network = (DM_Network*)dm->data; 783 784 PetscFunctionBegin; 785 ierr = PetscSectionGetOffset(network->DataSection,p,&offset);CHKERRQ(ierr); 786 *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata; 787 PetscFunctionReturn(0); 788 } 789 790 /*@ 791 DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector. 792 793 Not Collective 794 795 Input Parameters: 796 + dm - The DMNetwork object 797 - p - the edge/vertex point 798 799 Output Parameters: 800 . offset - the offset 801 802 Level: beginner 803 804 .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector 805 @*/ 806 PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset) 807 { 808 PetscErrorCode ierr; 809 DM_Network *network = (DM_Network*)dm->data; 810 811 PetscFunctionBegin; 812 ierr = PetscSectionGetOffset(network->plex->localSection,p,offset);CHKERRQ(ierr); 813 PetscFunctionReturn(0); 814 } 815 816 /*@ 817 DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector. 818 819 Not Collective 820 821 Input Parameters: 822 + dm - The DMNetwork object 823 - p - the edge/vertex point 824 825 Output Parameters: 826 . offsetg - the offset 827 828 Level: beginner 829 830 .seealso: DMNetworkGetVariableOffset, DMGetLocalVector 831 @*/ 832 PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg) 833 { 834 PetscErrorCode ierr; 835 DM_Network *network = (DM_Network*)dm->data; 836 837 PetscFunctionBegin; 838 ierr = PetscSectionGetOffset(network->plex->globalSection,p,offsetg);CHKERRQ(ierr); 839 if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */ 840 PetscFunctionReturn(0); 841 } 842 843 /*@ 844 DMNetworkGetComponentVariableOffset - Get the offset for accessing the variable associated with a component for the given vertex/edge from the local vector. 845 846 Not Collective 847 848 Input Parameters: 849 + dm - The DMNetwork object 850 . compnum - component number 851 - p - the edge/vertex point 852 853 Output Parameters: 854 . offset - the offset 855 856 Level: intermediate 857 858 .seealso: DMNetworkGetVariableGlobalOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables() 859 @*/ 860 PetscErrorCode DMNetworkGetComponentVariableOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset) 861 { 862 PetscErrorCode ierr; 863 DM_Network *network = (DM_Network*)dm->data; 864 PetscInt offsetp,offsetd; 865 DMNetworkComponentHeader header; 866 867 PetscFunctionBegin; 868 ierr = DMNetworkGetVariableOffset(dm,p,&offsetp);CHKERRQ(ierr); 869 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetd);CHKERRQ(ierr); 870 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd); 871 *offset = offsetp + header->offsetvarrel[compnum]; 872 PetscFunctionReturn(0); 873 } 874 875 /*@ 876 DMNetworkGetComponentVariableGlobalOffset - Get the global offset for accessing the variable associated with a component for the given vertex/edge from the local vector. 877 878 Not Collective 879 880 Input Parameters: 881 + dm - The DMNetwork object 882 . compnum - component number 883 - p - the edge/vertex point 884 885 Output Parameters: 886 . offsetg - the global offset 887 888 Level: intermediate 889 890 .seealso: DMNetworkGetVariableGlobalOffset(), DMNetworkGetComponentVariableOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables() 891 @*/ 892 PetscErrorCode DMNetworkGetComponentVariableGlobalOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg) 893 { 894 PetscErrorCode ierr; 895 DM_Network *network = (DM_Network*)dm->data; 896 PetscInt offsetp,offsetd; 897 DMNetworkComponentHeader header; 898 899 PetscFunctionBegin; 900 ierr = DMNetworkGetVariableGlobalOffset(dm,p,&offsetp);CHKERRQ(ierr); 901 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetd);CHKERRQ(ierr); 902 header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd); 903 *offsetg = offsetp + header->offsetvarrel[compnum]; 904 PetscFunctionReturn(0); 905 } 906 907 /*@ 908 DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector. 909 910 Not Collective 911 912 Input Parameters: 913 + dm - The DMNetwork object 914 - p - the edge point 915 916 Output Parameters: 917 . offset - the offset 918 919 Level: intermediate 920 921 .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector 922 @*/ 923 PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset) 924 { 925 PetscErrorCode ierr; 926 DM_Network *network = (DM_Network*)dm->data; 927 928 PetscFunctionBegin; 929 930 ierr = PetscSectionGetOffset(network->edge.DofSection,p,offset);CHKERRQ(ierr); 931 PetscFunctionReturn(0); 932 } 933 934 /*@ 935 DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector. 936 937 Not Collective 938 939 Input Parameters: 940 + dm - The DMNetwork object 941 - p - the vertex point 942 943 Output Parameters: 944 . offset - the offset 945 946 Level: intermediate 947 948 .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector 949 @*/ 950 PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset) 951 { 952 PetscErrorCode ierr; 953 DM_Network *network = (DM_Network*)dm->data; 954 955 PetscFunctionBegin; 956 957 p -= network->vStart; 958 959 ierr = PetscSectionGetOffset(network->vertex.DofSection,p,offset);CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 /*@ 963 DMNetworkAddNumVariables - Add number of variables associated with a given point. 964 965 Not Collective 966 967 Input Parameters: 968 + dm - The DMNetworkObject 969 . p - the vertex/edge point 970 - nvar - number of additional variables 971 972 Level: beginner 973 974 .seealso: DMNetworkSetNumVariables 975 @*/ 976 PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar) 977 { 978 PetscErrorCode ierr; 979 DM_Network *network = (DM_Network*)dm->data; 980 981 PetscFunctionBegin; 982 ierr = PetscSectionAddDof(network->DofSection,p,nvar);CHKERRQ(ierr); 983 PetscFunctionReturn(0); 984 } 985 986 /*@ 987 DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point. 988 989 Not Collective 990 991 Input Parameters: 992 + dm - The DMNetworkObject 993 - p - the vertex/edge point 994 995 Output Parameters: 996 . nvar - number of variables 997 998 Level: beginner 999 1000 .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables 1001 @*/ 1002 PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar) 1003 { 1004 PetscErrorCode ierr; 1005 DM_Network *network = (DM_Network*)dm->data; 1006 1007 PetscFunctionBegin; 1008 ierr = PetscSectionGetDof(network->DofSection,p,nvar);CHKERRQ(ierr); 1009 PetscFunctionReturn(0); 1010 } 1011 1012 /*@ 1013 DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point. 1014 1015 Not Collective 1016 1017 Input Parameters: 1018 + dm - The DMNetworkObject 1019 . p - the vertex/edge point 1020 - nvar - number of variables 1021 1022 Level: beginner 1023 1024 .seealso: DMNetworkAddNumVariables 1025 @*/ 1026 PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar) 1027 { 1028 PetscErrorCode ierr; 1029 DM_Network *network = (DM_Network*)dm->data; 1030 1031 PetscFunctionBegin; 1032 ierr = PetscSectionSetDof(network->DofSection,p,nvar);CHKERRQ(ierr); 1033 PetscFunctionReturn(0); 1034 } 1035 1036 /* Sets up the array that holds the data for all components and its associated section. This 1037 function is called during DMSetUp() */ 1038 PetscErrorCode DMNetworkComponentSetUp(DM dm) 1039 { 1040 PetscErrorCode ierr; 1041 DM_Network *network = (DM_Network*)dm->data; 1042 PetscInt arr_size,p,offset,offsetp,ncomp,i; 1043 DMNetworkComponentHeader header; 1044 DMNetworkComponentValue cvalue; 1045 DMNetworkComponentGenericDataType *componentdataarray; 1046 1047 PetscFunctionBegin; 1048 ierr = PetscSectionSetUp(network->DataSection);CHKERRQ(ierr); 1049 ierr = PetscSectionGetStorageSize(network->DataSection,&arr_size);CHKERRQ(ierr); 1050 ierr = PetscMalloc1(arr_size,&network->componentdataarray);CHKERRQ(ierr); 1051 componentdataarray = network->componentdataarray; 1052 for (p = network->pStart; p < network->pEnd; p++) { 1053 ierr = PetscSectionGetOffset(network->DataSection,p,&offsetp);CHKERRQ(ierr); 1054 /* Copy header */ 1055 header = &network->header[p]; 1056 ierr = PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));CHKERRQ(ierr); 1057 /* Copy data */ 1058 cvalue = &network->cvalue[p]; 1059 ncomp = header->ndata; 1060 for (i = 0; i < ncomp; i++) { 1061 offset = offsetp + network->dataheadersize + header->offset[i]; 1062 ierr = PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));CHKERRQ(ierr); 1063 } 1064 } 1065 PetscFunctionReturn(0); 1066 } 1067 1068 /* Sets up the section for dofs. This routine is called during DMSetUp() */ 1069 PetscErrorCode DMNetworkVariablesSetUp(DM dm) 1070 { 1071 PetscErrorCode ierr; 1072 DM_Network *network = (DM_Network*)dm->data; 1073 1074 PetscFunctionBegin; 1075 ierr = PetscSectionSetUp(network->DofSection);CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 /* 1080 DMNetworkGetComponentDataArray - Returns the component data array 1081 1082 Not Collective 1083 1084 Input Parameters: 1085 . dm - The DMNetwork Object 1086 1087 Output Parameters: 1088 . componentdataarray - array that holds data for all components 1089 1090 Level: intermediate 1091 1092 .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents 1093 */ 1094 PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray) 1095 { 1096 DM_Network *network = (DM_Network*)dm->data; 1097 1098 PetscFunctionBegin; 1099 *componentdataarray = network->componentdataarray; 1100 PetscFunctionReturn(0); 1101 } 1102 1103 /* Get a subsection from a range of points */ 1104 PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection) 1105 { 1106 PetscErrorCode ierr; 1107 PetscInt i, nvar; 1108 1109 PetscFunctionBegin; 1110 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);CHKERRQ(ierr); 1111 ierr = PetscSectionSetChart(*subsection, 0, pend - pstart);CHKERRQ(ierr); 1112 for (i = pstart; i < pend; i++) { 1113 ierr = PetscSectionGetDof(master,i,&nvar);CHKERRQ(ierr); 1114 ierr = PetscSectionSetDof(*subsection, i - pstart, nvar);CHKERRQ(ierr); 1115 } 1116 1117 ierr = PetscSectionSetUp(*subsection);CHKERRQ(ierr); 1118 PetscFunctionReturn(0); 1119 } 1120 1121 /* Create a submap of points with a GlobalToLocal structure */ 1122 PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map) 1123 { 1124 PetscErrorCode ierr; 1125 PetscInt i, *subpoints; 1126 1127 PetscFunctionBegin; 1128 /* Create index sets to map from "points" to "subpoints" */ 1129 ierr = PetscMalloc1(pend - pstart, &subpoints);CHKERRQ(ierr); 1130 for (i = pstart; i < pend; i++) { 1131 subpoints[i - pstart] = i; 1132 } 1133 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);CHKERRQ(ierr); 1134 ierr = PetscFree(subpoints);CHKERRQ(ierr); 1135 PetscFunctionReturn(0); 1136 } 1137 1138 /*@ 1139 DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute. 1140 1141 Collective 1142 1143 Input Parameters: 1144 . dm - The DMNetworkObject 1145 1146 Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are: 1147 1148 points = [0 1 2 3 4 5 6] 1149 1150 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]). 1151 1152 With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset. 1153 1154 Level: intermediate 1155 1156 @*/ 1157 PetscErrorCode DMNetworkAssembleGraphStructures(DM dm) 1158 { 1159 PetscErrorCode ierr; 1160 MPI_Comm comm; 1161 PetscMPIInt rank, size; 1162 DM_Network *network = (DM_Network*)dm->data; 1163 1164 PetscFunctionBegin; 1165 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1166 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1167 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1168 1169 /* Create maps for vertices and edges */ 1170 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 1171 ierr = DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);CHKERRQ(ierr); 1172 1173 /* Create local sub-sections */ 1174 ierr = DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);CHKERRQ(ierr); 1175 ierr = DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);CHKERRQ(ierr); 1176 1177 if (size > 1) { 1178 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 1179 1180 ierr = PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);CHKERRQ(ierr); 1181 ierr = PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);CHKERRQ(ierr); 1182 ierr = PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);CHKERRQ(ierr); 1183 } else { 1184 /* create structures for vertex */ 1185 ierr = PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);CHKERRQ(ierr); 1186 /* create structures for edge */ 1187 ierr = PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);CHKERRQ(ierr); 1188 } 1189 1190 /* Add viewers */ 1191 ierr = PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");CHKERRQ(ierr); 1192 ierr = PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");CHKERRQ(ierr); 1193 ierr = PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");CHKERRQ(ierr); 1194 ierr = PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /*@ 1199 DMNetworkDistribute - Distributes the network and moves associated component data. 1200 1201 Collective 1202 1203 Input Parameter: 1204 + DM - the DMNetwork object 1205 - overlap - The overlap of partitions, 0 is the default 1206 1207 Notes: 1208 Distributes the network with <overlap>-overlapping partitioning of the edges. 1209 1210 Level: intermediate 1211 1212 .seealso: DMNetworkCreate 1213 @*/ 1214 PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap) 1215 { 1216 MPI_Comm comm; 1217 PetscErrorCode ierr; 1218 PetscMPIInt size; 1219 DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data); 1220 DM_Network *newDMnetwork; 1221 PetscSF pointsf=NULL; 1222 DM newDM; 1223 PetscInt j,e,v,offset,*subnetvtx; 1224 PetscPartitioner part; 1225 DMNetworkComponentHeader header; 1226 1227 PetscFunctionBegin; 1228 ierr = PetscObjectGetComm((PetscObject)*dm,&comm);CHKERRQ(ierr); 1229 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1230 if (size == 1) PetscFunctionReturn(0); 1231 1232 ierr = DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);CHKERRQ(ierr); 1233 newDMnetwork = (DM_Network*)newDM->data; 1234 newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType); 1235 1236 /* Enable runtime options for petscpartitioner */ 1237 ierr = DMPlexGetPartitioner(oldDMnetwork->plex,&part);CHKERRQ(ierr); 1238 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 1239 1240 /* Distribute plex dm and dof section */ 1241 ierr = DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);CHKERRQ(ierr); 1242 1243 /* Distribute dof section */ 1244 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);CHKERRQ(ierr); 1245 ierr = PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);CHKERRQ(ierr); 1246 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);CHKERRQ(ierr); 1247 1248 /* Distribute data and associated section */ 1249 ierr = DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);CHKERRQ(ierr); 1250 1251 ierr = PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);CHKERRQ(ierr); 1252 ierr = DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);CHKERRQ(ierr); 1253 ierr = DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);CHKERRQ(ierr); 1254 newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart; 1255 newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart; 1256 newDMnetwork->NVertices = oldDMnetwork->NVertices; 1257 newDMnetwork->NEdges = oldDMnetwork->NEdges; 1258 1259 /* Set Dof section as the section for dm */ 1260 ierr = DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);CHKERRQ(ierr); 1261 ierr = DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);CHKERRQ(ierr); 1262 1263 /* Set up subnetwork info in the newDM */ 1264 newDMnetwork->nsubnet = oldDMnetwork->nsubnet; 1265 newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet; 1266 ierr = PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);CHKERRQ(ierr); 1267 /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already 1268 calculated in DMNetworkLayoutSetUp() 1269 */ 1270 for(j=0; j < newDMnetwork->nsubnet; j++) { 1271 newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx; 1272 newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge; 1273 } 1274 1275 for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) { 1276 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);CHKERRQ(ierr); 1277 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1278 newDMnetwork->subnet[header->subnetid].nedge++; 1279 } 1280 1281 for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) { 1282 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);CHKERRQ(ierr); 1283 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1284 newDMnetwork->subnet[header->subnetid].nvtx++; 1285 } 1286 1287 /* Now create the vertices and edge arrays for the subnetworks */ 1288 ierr = PetscCalloc1(newDMnetwork->vEnd-newDMnetwork->vStart,&newDMnetwork->subnetvtx);CHKERRQ(ierr); 1289 subnetvtx = newDMnetwork->subnetvtx; 1290 1291 for (j=0; j<newDMnetwork->nsubnet; j++) { 1292 ierr = PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);CHKERRQ(ierr); 1293 newDMnetwork->subnet[j].vertices = subnetvtx; 1294 subnetvtx += newDMnetwork->subnet[j].nvtx; 1295 1296 /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop. 1297 These get updated when the vertices and edges are added. */ 1298 newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0; 1299 } 1300 1301 /* Set the vertices and edges in each subnetwork */ 1302 for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) { 1303 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);CHKERRQ(ierr); 1304 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1305 newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e; 1306 } 1307 1308 for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) { 1309 ierr = PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);CHKERRQ(ierr); 1310 header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);CHKERRQ(ierr); 1311 newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v; 1312 } 1313 1314 newDM->setupcalled = (*dm)->setupcalled; 1315 newDMnetwork->distributecalled = PETSC_TRUE; 1316 1317 /* Destroy point SF */ 1318 ierr = PetscSFDestroy(&pointsf);CHKERRQ(ierr); 1319 1320 ierr = DMDestroy(dm);CHKERRQ(ierr); 1321 *dm = newDM; 1322 PetscFunctionReturn(0); 1323 } 1324 1325 /*@C 1326 PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering. 1327 1328 Input Parameters: 1329 + masterSF - the original SF structure 1330 - map - a ISLocalToGlobal mapping that contains the subset of points 1331 1332 Output Parameters: 1333 . subSF - a subset of the masterSF for the desired subset. 1334 */ 1335 1336 PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) { 1337 1338 PetscErrorCode ierr; 1339 PetscInt nroots, nleaves, *ilocal_sub; 1340 PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0; 1341 PetscInt *local_points, *remote_points; 1342 PetscSFNode *iremote_sub; 1343 const PetscInt *ilocal; 1344 const PetscSFNode *iremote; 1345 1346 PetscFunctionBegin; 1347 ierr = PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 1348 1349 /* Look for leaves that pertain to the subset of points. Get the local ordering */ 1350 ierr = PetscMalloc1(nleaves,&ilocal_map);CHKERRQ(ierr); 1351 ierr = ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);CHKERRQ(ierr); 1352 for (i = 0; i < nleaves; i++) { 1353 if (ilocal_map[i] != -1) nleaves_sub += 1; 1354 } 1355 /* Re-number ilocal with subset numbering. Need information from roots */ 1356 ierr = PetscMalloc2(nroots,&local_points,nroots,&remote_points);CHKERRQ(ierr); 1357 for (i = 0; i < nroots; i++) local_points[i] = i; 1358 ierr = ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);CHKERRQ(ierr); 1359 ierr = PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);CHKERRQ(ierr); 1360 ierr = PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);CHKERRQ(ierr); 1361 /* Fill up graph using local (that is, local to the subset) numbering. */ 1362 ierr = PetscMalloc1(nleaves_sub,&ilocal_sub);CHKERRQ(ierr); 1363 ierr = PetscMalloc1(nleaves_sub,&iremote_sub);CHKERRQ(ierr); 1364 nleaves_sub = 0; 1365 for (i = 0; i < nleaves; i++) { 1366 if (ilocal_map[i] != -1) { 1367 ilocal_sub[nleaves_sub] = ilocal_map[i]; 1368 iremote_sub[nleaves_sub].rank = iremote[i].rank; 1369 iremote_sub[nleaves_sub].index = remote_points[ilocal[i]]; 1370 nleaves_sub += 1; 1371 } 1372 } 1373 ierr = PetscFree2(local_points,remote_points);CHKERRQ(ierr); 1374 ierr = ISLocalToGlobalMappingGetSize(map,&nroots_sub);CHKERRQ(ierr); 1375 1376 /* Create new subSF */ 1377 ierr = PetscSFCreate(PETSC_COMM_WORLD,subSF);CHKERRQ(ierr); 1378 ierr = PetscSFSetFromOptions(*subSF);CHKERRQ(ierr); 1379 ierr = PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);CHKERRQ(ierr); 1380 ierr = PetscFree(ilocal_map);CHKERRQ(ierr); 1381 ierr = PetscFree(iremote_sub);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 /*@C 1386 DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point 1387 1388 Not Collective 1389 1390 Input Parameters: 1391 + dm - The DMNetwork object 1392 - p - the vertex point 1393 1394 Output Parameters: 1395 + nedges - number of edges connected to this vertex point 1396 - edges - List of edge points 1397 1398 Level: beginner 1399 1400 Fortran Notes: 1401 Since it returns an array, this routine is only available in Fortran 90, and you must 1402 include petsc.h90 in your code. 1403 1404 .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices 1405 @*/ 1406 PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[]) 1407 { 1408 PetscErrorCode ierr; 1409 DM_Network *network = (DM_Network*)dm->data; 1410 1411 PetscFunctionBegin; 1412 ierr = DMPlexGetSupportSize(network->plex,vertex,nedges);CHKERRQ(ierr); 1413 ierr = DMPlexGetSupport(network->plex,vertex,edges);CHKERRQ(ierr); 1414 PetscFunctionReturn(0); 1415 } 1416 1417 /*@C 1418 DMNetworkGetConnectedVertices - Return the connected vertices for this edge point 1419 1420 Not Collective 1421 1422 Input Parameters: 1423 + dm - The DMNetwork object 1424 - p - the edge point 1425 1426 Output Parameters: 1427 . vertices - vertices connected to this edge 1428 1429 Level: beginner 1430 1431 Fortran Notes: 1432 Since it returns an array, this routine is only available in Fortran 90, and you must 1433 include petsc.h90 in your code. 1434 1435 .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges 1436 @*/ 1437 PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[]) 1438 { 1439 PetscErrorCode ierr; 1440 DM_Network *network = (DM_Network*)dm->data; 1441 1442 PetscFunctionBegin; 1443 ierr = DMPlexGetCone(network->plex,edge,vertices);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 /*@ 1448 DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex 1449 1450 Not Collective 1451 1452 Input Parameters: 1453 + dm - The DMNetwork object 1454 - p - the vertex point 1455 1456 Output Parameter: 1457 . isghost - TRUE if the vertex is a ghost point 1458 1459 Level: beginner 1460 1461 .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange 1462 @*/ 1463 PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost) 1464 { 1465 PetscErrorCode ierr; 1466 DM_Network *network = (DM_Network*)dm->data; 1467 PetscInt offsetg; 1468 PetscSection sectiong; 1469 1470 PetscFunctionBegin; 1471 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 1472 *isghost = PETSC_FALSE; 1473 ierr = DMGetGlobalSection(network->plex,§iong);CHKERRQ(ierr); 1474 ierr = PetscSectionGetOffset(sectiong,p,&offsetg);CHKERRQ(ierr); 1475 if (offsetg < 0) *isghost = PETSC_TRUE; 1476 PetscFunctionReturn(0); 1477 } 1478 1479 PetscErrorCode DMSetUp_Network(DM dm) 1480 { 1481 PetscErrorCode ierr; 1482 DM_Network *network=(DM_Network*)dm->data; 1483 1484 PetscFunctionBegin; 1485 ierr = DMNetworkComponentSetUp(dm);CHKERRQ(ierr); 1486 ierr = DMNetworkVariablesSetUp(dm);CHKERRQ(ierr); 1487 1488 ierr = DMSetLocalSection(network->plex,network->DofSection);CHKERRQ(ierr); 1489 ierr = DMGetGlobalSection(network->plex,&network->GlobalDofSection);CHKERRQ(ierr); 1490 1491 dm->setupcalled = PETSC_TRUE; 1492 ierr = DMViewFromOptions(dm,NULL,"-dm_view");CHKERRQ(ierr); 1493 PetscFunctionReturn(0); 1494 } 1495 1496 /*@ 1497 DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices 1498 -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)? 1499 1500 Collective 1501 1502 Input Parameters: 1503 + dm - The DMNetwork object 1504 . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges 1505 - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices 1506 1507 Level: intermediate 1508 1509 @*/ 1510 PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg) 1511 { 1512 DM_Network *network=(DM_Network*)dm->data; 1513 PetscErrorCode ierr; 1514 PetscInt nVertices = network->nVertices; 1515 1516 PetscFunctionBegin; 1517 network->userEdgeJacobian = eflg; 1518 network->userVertexJacobian = vflg; 1519 1520 if (eflg && !network->Je) { 1521 ierr = PetscCalloc1(3*network->nEdges,&network->Je);CHKERRQ(ierr); 1522 } 1523 1524 if (vflg && !network->Jv && nVertices) { 1525 PetscInt i,*vptr,nedges,vStart=network->vStart; 1526 PetscInt nedges_total; 1527 const PetscInt *edges; 1528 1529 /* count nvertex_total */ 1530 nedges_total = 0; 1531 ierr = PetscMalloc1(nVertices+1,&vptr);CHKERRQ(ierr); 1532 1533 vptr[0] = 0; 1534 for (i=0; i<nVertices; i++) { 1535 ierr = DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);CHKERRQ(ierr); 1536 nedges_total += nedges; 1537 vptr[i+1] = vptr[i] + 2*nedges + 1; 1538 } 1539 1540 ierr = PetscCalloc1(2*nedges_total+nVertices,&network->Jv);CHKERRQ(ierr); 1541 network->Jvptr = vptr; 1542 } 1543 PetscFunctionReturn(0); 1544 } 1545 1546 /*@ 1547 DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network 1548 1549 Not Collective 1550 1551 Input Parameters: 1552 + dm - The DMNetwork object 1553 . p - the edge point 1554 - J - array (size = 3) of Jacobian submatrices for this edge point: 1555 J[0]: this edge 1556 J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices() 1557 1558 Level: advanced 1559 1560 .seealso: DMNetworkVertexSetMatrix 1561 @*/ 1562 PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[]) 1563 { 1564 DM_Network *network=(DM_Network*)dm->data; 1565 1566 PetscFunctionBegin; 1567 if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix"); 1568 1569 if (J) { 1570 network->Je[3*p] = J[0]; 1571 network->Je[3*p+1] = J[1]; 1572 network->Je[3*p+2] = J[2]; 1573 } 1574 PetscFunctionReturn(0); 1575 } 1576 1577 /*@ 1578 DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network 1579 1580 Not Collective 1581 1582 Input Parameters: 1583 + dm - The DMNetwork object 1584 . p - the vertex point 1585 - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point: 1586 J[0]: this vertex 1587 J[1+2*i]: i-th supporting edge 1588 J[1+2*i+1]: i-th connected vertex 1589 1590 Level: advanced 1591 1592 .seealso: DMNetworkEdgeSetMatrix 1593 @*/ 1594 PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[]) 1595 { 1596 PetscErrorCode ierr; 1597 DM_Network *network=(DM_Network*)dm->data; 1598 PetscInt i,*vptr,nedges,vStart=network->vStart; 1599 const PetscInt *edges; 1600 1601 PetscFunctionBegin; 1602 if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix"); 1603 1604 if (J) { 1605 vptr = network->Jvptr; 1606 network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */ 1607 1608 /* Set Jacobian for each supporting edge and connected vertex */ 1609 ierr = DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);CHKERRQ(ierr); 1610 for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i]; 1611 } 1612 PetscFunctionReturn(0); 1613 } 1614 1615 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 1616 { 1617 PetscErrorCode ierr; 1618 PetscInt j; 1619 PetscScalar val=(PetscScalar)ncols; 1620 1621 PetscFunctionBegin; 1622 if (!ghost) { 1623 for (j=0; j<nrows; j++) { 1624 ierr = VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1625 } 1626 } else { 1627 for (j=0; j<nrows; j++) { 1628 ierr = VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1629 } 1630 } 1631 PetscFunctionReturn(0); 1632 } 1633 1634 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 1635 { 1636 PetscErrorCode ierr; 1637 PetscInt j,ncols_u; 1638 PetscScalar val; 1639 1640 PetscFunctionBegin; 1641 if (!ghost) { 1642 for (j=0; j<nrows; j++) { 1643 ierr = MatGetRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1644 val = (PetscScalar)ncols_u; 1645 ierr = VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1646 ierr = MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1647 } 1648 } else { 1649 for (j=0; j<nrows; j++) { 1650 ierr = MatGetRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1651 val = (PetscScalar)ncols_u; 1652 ierr = VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);CHKERRQ(ierr); 1653 ierr = MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);CHKERRQ(ierr); 1654 } 1655 } 1656 PetscFunctionReturn(0); 1657 } 1658 1659 PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz) 1660 { 1661 PetscErrorCode ierr; 1662 1663 PetscFunctionBegin; 1664 if (Ju) { 1665 ierr = MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);CHKERRQ(ierr); 1666 } else { 1667 ierr = MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);CHKERRQ(ierr); 1668 } 1669 PetscFunctionReturn(0); 1670 } 1671 1672 PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 1673 { 1674 PetscErrorCode ierr; 1675 PetscInt j,*cols; 1676 PetscScalar *zeros; 1677 1678 PetscFunctionBegin; 1679 ierr = PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);CHKERRQ(ierr); 1680 for (j=0; j<ncols; j++) cols[j] = j+ cstart; 1681 ierr = MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);CHKERRQ(ierr); 1682 ierr = PetscFree2(cols,zeros);CHKERRQ(ierr); 1683 PetscFunctionReturn(0); 1684 } 1685 1686 PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 1687 { 1688 PetscErrorCode ierr; 1689 PetscInt j,M,N,row,col,ncols_u; 1690 const PetscInt *cols; 1691 PetscScalar zero=0.0; 1692 1693 PetscFunctionBegin; 1694 ierr = MatGetSize(Ju,&M,&N);CHKERRQ(ierr); 1695 if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N); 1696 1697 for (row=0; row<nrows; row++) { 1698 ierr = MatGetRow(Ju,row,&ncols_u,&cols,NULL);CHKERRQ(ierr); 1699 for (j=0; j<ncols_u; j++) { 1700 col = cols[j] + cstart; 1701 ierr = MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);CHKERRQ(ierr); 1702 } 1703 ierr = MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);CHKERRQ(ierr); 1704 } 1705 PetscFunctionReturn(0); 1706 } 1707 1708 PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J) 1709 { 1710 PetscErrorCode ierr; 1711 1712 PetscFunctionBegin; 1713 if (Ju) { 1714 ierr = MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 1715 } else { 1716 ierr = MatSetDenseblock_private(nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 1717 } 1718 PetscFunctionReturn(0); 1719 } 1720 1721 /* 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. 1722 */ 1723 PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog) 1724 { 1725 PetscErrorCode ierr; 1726 PetscInt i,size,dof; 1727 PetscInt *glob2loc; 1728 1729 PetscFunctionBegin; 1730 ierr = PetscSectionGetStorageSize(localsec,&size);CHKERRQ(ierr); 1731 ierr = PetscMalloc1(size,&glob2loc);CHKERRQ(ierr); 1732 1733 for (i = 0; i < size; i++) { 1734 ierr = PetscSectionGetOffset(globalsec,i,&dof);CHKERRQ(ierr); 1735 dof = (dof >= 0) ? dof : -(dof + 1); 1736 glob2loc[i] = dof; 1737 } 1738 1739 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1740 #if 0 1741 ierr = PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1742 #endif 1743 PetscFunctionReturn(0); 1744 } 1745 1746 #include <petsc/private/matimpl.h> 1747 1748 PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J) 1749 { 1750 PetscErrorCode ierr; 1751 DM_Network *network = (DM_Network*)dm->data; 1752 PetscMPIInt rank, size; 1753 PetscInt eDof,vDof; 1754 Mat j11,j12,j21,j22,bA[2][2]; 1755 MPI_Comm comm; 1756 ISLocalToGlobalMapping eISMap,vISMap; 1757 1758 PetscFunctionBegin; 1759 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1760 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1761 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1762 1763 ierr = PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);CHKERRQ(ierr); 1764 ierr = PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);CHKERRQ(ierr); 1765 1766 ierr = MatCreate(comm, &j11);CHKERRQ(ierr); 1767 ierr = MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1768 ierr = MatSetType(j11, MATMPIAIJ);CHKERRQ(ierr); 1769 1770 ierr = MatCreate(comm, &j12);CHKERRQ(ierr); 1771 ierr = MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);CHKERRQ(ierr); 1772 ierr = MatSetType(j12, MATMPIAIJ);CHKERRQ(ierr); 1773 1774 ierr = MatCreate(comm, &j21);CHKERRQ(ierr); 1775 ierr = MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1776 ierr = MatSetType(j21, MATMPIAIJ);CHKERRQ(ierr); 1777 1778 ierr = MatCreate(comm, &j22);CHKERRQ(ierr); 1779 ierr = MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1780 ierr = MatSetType(j22, MATMPIAIJ);CHKERRQ(ierr); 1781 1782 bA[0][0] = j11; 1783 bA[0][1] = j12; 1784 bA[1][0] = j21; 1785 bA[1][1] = j22; 1786 1787 ierr = CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);CHKERRQ(ierr); 1788 ierr = CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);CHKERRQ(ierr); 1789 1790 ierr = MatSetLocalToGlobalMapping(j11,eISMap,eISMap);CHKERRQ(ierr); 1791 ierr = MatSetLocalToGlobalMapping(j12,eISMap,vISMap);CHKERRQ(ierr); 1792 ierr = MatSetLocalToGlobalMapping(j21,vISMap,eISMap);CHKERRQ(ierr); 1793 ierr = MatSetLocalToGlobalMapping(j22,vISMap,vISMap);CHKERRQ(ierr); 1794 1795 ierr = MatSetUp(j11);CHKERRQ(ierr); 1796 ierr = MatSetUp(j12);CHKERRQ(ierr); 1797 ierr = MatSetUp(j21);CHKERRQ(ierr); 1798 ierr = MatSetUp(j22);CHKERRQ(ierr); 1799 1800 ierr = MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);CHKERRQ(ierr); 1801 ierr = MatSetUp(*J);CHKERRQ(ierr); 1802 ierr = MatNestSetVecType(*J,VECNEST);CHKERRQ(ierr); 1803 ierr = MatDestroy(&j11);CHKERRQ(ierr); 1804 ierr = MatDestroy(&j12);CHKERRQ(ierr); 1805 ierr = MatDestroy(&j21);CHKERRQ(ierr); 1806 ierr = MatDestroy(&j22);CHKERRQ(ierr); 1807 1808 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1809 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1810 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1811 1812 /* Free structures */ 1813 ierr = ISLocalToGlobalMappingDestroy(&eISMap);CHKERRQ(ierr); 1814 ierr = ISLocalToGlobalMappingDestroy(&vISMap);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J) 1819 { 1820 PetscErrorCode ierr; 1821 DM_Network *network = (DM_Network*)dm->data; 1822 PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize; 1823 PetscInt cstart,ncols,j,e,v; 1824 PetscBool ghost,ghost_vc,ghost2,isNest; 1825 Mat Juser; 1826 PetscSection sectionGlobal; 1827 PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */ 1828 const PetscInt *edges,*cone; 1829 MPI_Comm comm; 1830 MatType mtype; 1831 Vec vd_nz,vo_nz; 1832 PetscInt *dnnz,*onnz; 1833 PetscScalar *vdnz,*vonz; 1834 1835 PetscFunctionBegin; 1836 mtype = dm->mattype; 1837 ierr = PetscStrcmp(mtype,MATNEST,&isNest);CHKERRQ(ierr); 1838 if (isNest) { 1839 ierr = DMCreateMatrix_Network_Nest(dm,J);CHKERRQ(ierr); 1840 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 1841 PetscFunctionReturn(0); 1842 } 1843 1844 if (!network->userEdgeJacobian && !network->userVertexJacobian) { 1845 /* user does not provide Jacobian blocks */ 1846 ierr = DMCreateMatrix_Plex(network->plex,J);CHKERRQ(ierr); 1847 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 1848 PetscFunctionReturn(0); 1849 } 1850 1851 ierr = MatCreate(PetscObjectComm((PetscObject)dm),J);CHKERRQ(ierr); 1852 ierr = DMGetGlobalSection(network->plex,§ionGlobal);CHKERRQ(ierr); 1853 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);CHKERRQ(ierr); 1854 ierr = MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1855 1856 ierr = MatSetType(*J,MATAIJ);CHKERRQ(ierr); 1857 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 1858 1859 /* (1) Set matrix preallocation */ 1860 /*------------------------------*/ 1861 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1862 ierr = VecCreate(comm,&vd_nz);CHKERRQ(ierr); 1863 ierr = VecSetSizes(vd_nz,localSize,PETSC_DECIDE);CHKERRQ(ierr); 1864 ierr = VecSetFromOptions(vd_nz);CHKERRQ(ierr); 1865 ierr = VecSet(vd_nz,0.0);CHKERRQ(ierr); 1866 ierr = VecDuplicate(vd_nz,&vo_nz);CHKERRQ(ierr); 1867 1868 /* Set preallocation for edges */ 1869 /*-----------------------------*/ 1870 ierr = DMNetworkGetEdgeRange(dm,&eStart,&eEnd);CHKERRQ(ierr); 1871 1872 ierr = PetscMalloc1(localSize,&rows);CHKERRQ(ierr); 1873 for (e=eStart; e<eEnd; e++) { 1874 /* Get row indices */ 1875 ierr = DMNetworkGetVariableGlobalOffset(dm,e,&rstart);CHKERRQ(ierr); 1876 ierr = DMNetworkGetNumVariables(dm,e,&nrows);CHKERRQ(ierr); 1877 if (nrows) { 1878 for (j=0; j<nrows; j++) rows[j] = j + rstart; 1879 1880 /* Set preallocation for conntected vertices */ 1881 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 1882 for (v=0; v<2; v++) { 1883 ierr = DMNetworkGetNumVariables(dm,cone[v],&ncols);CHKERRQ(ierr); 1884 1885 if (network->Je) { 1886 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 1887 } else Juser = NULL; 1888 ierr = DMNetworkIsGhostVertex(dm,cone[v],&ghost);CHKERRQ(ierr); 1889 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 1890 } 1891 1892 /* Set preallocation for edge self */ 1893 cstart = rstart; 1894 if (network->Je) { 1895 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 1896 } else Juser = NULL; 1897 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 1898 } 1899 } 1900 1901 /* Set preallocation for vertices */ 1902 /*--------------------------------*/ 1903 ierr = DMNetworkGetVertexRange(dm,&vStart,&vEnd);CHKERRQ(ierr); 1904 if (vEnd - vStart) vptr = network->Jvptr; 1905 1906 for (v=vStart; v<vEnd; v++) { 1907 /* Get row indices */ 1908 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&rstart);CHKERRQ(ierr); 1909 ierr = DMNetworkGetNumVariables(dm,v,&nrows);CHKERRQ(ierr); 1910 if (!nrows) continue; 1911 1912 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1913 if (ghost) { 1914 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 1915 } else { 1916 rows_v = rows; 1917 } 1918 1919 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 1920 1921 /* Get supporting edges and connected vertices */ 1922 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 1923 1924 for (e=0; e<nedges; e++) { 1925 /* Supporting edges */ 1926 ierr = DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);CHKERRQ(ierr); 1927 ierr = DMNetworkGetNumVariables(dm,edges[e],&ncols);CHKERRQ(ierr); 1928 1929 if (network->Jv) { 1930 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 1931 } else Juser = NULL; 1932 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 1933 1934 /* Connected vertices */ 1935 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 1936 vc = (v == cone[0]) ? cone[1]:cone[0]; 1937 ierr = DMNetworkIsGhostVertex(dm,vc,&ghost_vc);CHKERRQ(ierr); 1938 1939 ierr = DMNetworkGetNumVariables(dm,vc,&ncols);CHKERRQ(ierr); 1940 1941 if (network->Jv) { 1942 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 1943 } else Juser = NULL; 1944 if (ghost_vc||ghost) { 1945 ghost2 = PETSC_TRUE; 1946 } else { 1947 ghost2 = PETSC_FALSE; 1948 } 1949 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);CHKERRQ(ierr); 1950 } 1951 1952 /* Set preallocation for vertex self */ 1953 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1954 if (!ghost) { 1955 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&cstart);CHKERRQ(ierr); 1956 if (network->Jv) { 1957 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 1958 } else Juser = NULL; 1959 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 1960 } 1961 if (ghost) { 1962 ierr = PetscFree(rows_v);CHKERRQ(ierr); 1963 } 1964 } 1965 1966 ierr = VecAssemblyBegin(vd_nz);CHKERRQ(ierr); 1967 ierr = VecAssemblyBegin(vo_nz);CHKERRQ(ierr); 1968 1969 ierr = PetscMalloc2(localSize,&dnnz,localSize,&onnz);CHKERRQ(ierr); 1970 1971 ierr = VecAssemblyEnd(vd_nz);CHKERRQ(ierr); 1972 ierr = VecAssemblyEnd(vo_nz);CHKERRQ(ierr); 1973 1974 ierr = VecGetArray(vd_nz,&vdnz);CHKERRQ(ierr); 1975 ierr = VecGetArray(vo_nz,&vonz);CHKERRQ(ierr); 1976 for (j=0; j<localSize; j++) { 1977 dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]); 1978 onnz[j] = (PetscInt)PetscRealPart(vonz[j]); 1979 } 1980 ierr = VecRestoreArray(vd_nz,&vdnz);CHKERRQ(ierr); 1981 ierr = VecRestoreArray(vo_nz,&vonz);CHKERRQ(ierr); 1982 ierr = VecDestroy(&vd_nz);CHKERRQ(ierr); 1983 ierr = VecDestroy(&vo_nz);CHKERRQ(ierr); 1984 1985 ierr = MatSeqAIJSetPreallocation(*J,0,dnnz);CHKERRQ(ierr); 1986 ierr = MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);CHKERRQ(ierr); 1987 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1988 1989 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 1990 1991 /* (2) Set matrix entries for edges */ 1992 /*----------------------------------*/ 1993 for (e=eStart; e<eEnd; e++) { 1994 /* Get row indices */ 1995 ierr = DMNetworkGetVariableGlobalOffset(dm,e,&rstart);CHKERRQ(ierr); 1996 ierr = DMNetworkGetNumVariables(dm,e,&nrows);CHKERRQ(ierr); 1997 if (nrows) { 1998 for (j=0; j<nrows; j++) rows[j] = j + rstart; 1999 2000 /* Set matrix entries for conntected vertices */ 2001 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 2002 for (v=0; v<2; v++) { 2003 ierr = DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);CHKERRQ(ierr); 2004 ierr = DMNetworkGetNumVariables(dm,cone[v],&ncols);CHKERRQ(ierr); 2005 2006 if (network->Je) { 2007 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 2008 } else Juser = NULL; 2009 ierr = MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 2010 } 2011 2012 /* Set matrix entries for edge self */ 2013 cstart = rstart; 2014 if (network->Je) { 2015 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 2016 } else Juser = NULL; 2017 ierr = MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);CHKERRQ(ierr); 2018 } 2019 } 2020 2021 /* Set matrix entries for vertices */ 2022 /*---------------------------------*/ 2023 for (v=vStart; v<vEnd; v++) { 2024 /* Get row indices */ 2025 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&rstart);CHKERRQ(ierr); 2026 ierr = DMNetworkGetNumVariables(dm,v,&nrows);CHKERRQ(ierr); 2027 if (!nrows) continue; 2028 2029 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 2030 if (ghost) { 2031 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 2032 } else { 2033 rows_v = rows; 2034 } 2035 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 2036 2037 /* Get supporting edges and connected vertices */ 2038 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 2039 2040 for (e=0; e<nedges; e++) { 2041 /* Supporting edges */ 2042 ierr = DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);CHKERRQ(ierr); 2043 ierr = DMNetworkGetNumVariables(dm,edges[e],&ncols);CHKERRQ(ierr); 2044 2045 if (network->Jv) { 2046 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 2047 } else Juser = NULL; 2048 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 2049 2050 /* Connected vertices */ 2051 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 2052 vc = (v == cone[0]) ? cone[1]:cone[0]; 2053 2054 ierr = DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);CHKERRQ(ierr); 2055 ierr = DMNetworkGetNumVariables(dm,vc,&ncols);CHKERRQ(ierr); 2056 2057 if (network->Jv) { 2058 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 2059 } else Juser = NULL; 2060 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 2061 } 2062 2063 /* Set matrix entries for vertex self */ 2064 if (!ghost) { 2065 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&cstart);CHKERRQ(ierr); 2066 if (network->Jv) { 2067 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 2068 } else Juser = NULL; 2069 ierr = MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);CHKERRQ(ierr); 2070 } 2071 if (ghost) { 2072 ierr = PetscFree(rows_v);CHKERRQ(ierr); 2073 } 2074 } 2075 ierr = PetscFree(rows);CHKERRQ(ierr); 2076 2077 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2078 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2079 2080 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 2081 PetscFunctionReturn(0); 2082 } 2083 2084 PetscErrorCode DMDestroy_Network(DM dm) 2085 { 2086 PetscErrorCode ierr; 2087 DM_Network *network = (DM_Network*)dm->data; 2088 PetscInt j; 2089 2090 PetscFunctionBegin; 2091 if (--network->refct > 0) PetscFunctionReturn(0); 2092 if (network->Je) { 2093 ierr = PetscFree(network->Je);CHKERRQ(ierr); 2094 } 2095 if (network->Jv) { 2096 ierr = PetscFree(network->Jvptr);CHKERRQ(ierr); 2097 ierr = PetscFree(network->Jv);CHKERRQ(ierr); 2098 } 2099 2100 ierr = ISLocalToGlobalMappingDestroy(&network->vertex.mapping);CHKERRQ(ierr); 2101 ierr = PetscSectionDestroy(&network->vertex.DofSection);CHKERRQ(ierr); 2102 ierr = PetscSectionDestroy(&network->vertex.GlobalDofSection);CHKERRQ(ierr); 2103 if (network->vltog) { 2104 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2105 } 2106 if (network->vertex.sf) { 2107 ierr = PetscSFDestroy(&network->vertex.sf);CHKERRQ(ierr); 2108 } 2109 /* edge */ 2110 ierr = ISLocalToGlobalMappingDestroy(&network->edge.mapping);CHKERRQ(ierr); 2111 ierr = PetscSectionDestroy(&network->edge.DofSection);CHKERRQ(ierr); 2112 ierr = PetscSectionDestroy(&network->edge.GlobalDofSection);CHKERRQ(ierr); 2113 if (network->edge.sf) { 2114 ierr = PetscSFDestroy(&network->edge.sf);CHKERRQ(ierr); 2115 } 2116 ierr = DMDestroy(&network->plex);CHKERRQ(ierr); 2117 ierr = PetscSectionDestroy(&network->DataSection);CHKERRQ(ierr); 2118 ierr = PetscSectionDestroy(&network->DofSection);CHKERRQ(ierr); 2119 2120 for(j=0; j<network->nsubnet; j++) { 2121 ierr = PetscFree(network->subnet[j].edges);CHKERRQ(ierr); 2122 } 2123 ierr = PetscFree(network->subnetvtx);CHKERRQ(ierr); 2124 2125 ierr = PetscFree(network->subnet);CHKERRQ(ierr); 2126 ierr = PetscFree(network->componentdataarray);CHKERRQ(ierr); 2127 ierr = PetscFree2(network->header,network->cvalue);CHKERRQ(ierr); 2128 ierr = PetscFree(network);CHKERRQ(ierr); 2129 PetscFunctionReturn(0); 2130 } 2131 2132 PetscErrorCode DMView_Network(DM dm,PetscViewer viewer) 2133 { 2134 PetscErrorCode ierr; 2135 DM_Network *network = (DM_Network*)dm->data; 2136 PetscBool iascii; 2137 PetscMPIInt rank; 2138 PetscInt p,nsubnet; 2139 2140 PetscFunctionBegin; 2141 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 2142 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 2143 PetscValidHeaderSpecific(dm,DM_CLASSID, 1); 2144 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 2145 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2146 if (iascii) { 2147 const PetscInt *cone,*vtx,*edges; 2148 PetscInt vfrom,vto,i,j,nv,ne; 2149 2150 nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */ 2151 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 2152 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);CHKERRQ(ierr); 2153 2154 for (i=0; i<nsubnet; i++) { 2155 ierr = DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2156 if (ne) { 2157 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);CHKERRQ(ierr); 2158 for (j=0; j<ne; j++) { 2159 p = edges[j]; 2160 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2161 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2162 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2163 ierr = DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);CHKERRQ(ierr); 2164 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2165 } 2166 } 2167 } 2168 /* Coupling subnets */ 2169 nsubnet = network->nsubnet; 2170 for (; i<nsubnet; i++) { 2171 ierr = DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2172 if (ne) { 2173 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);CHKERRQ(ierr); 2174 for (j=0; j<ne; j++) { 2175 p = edges[j]; 2176 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2177 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2178 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2179 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2180 } 2181 } 2182 } 2183 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2184 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 2185 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name); 2186 PetscFunctionReturn(0); 2187 } 2188 2189 PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l) 2190 { 2191 PetscErrorCode ierr; 2192 DM_Network *network = (DM_Network*)dm->data; 2193 2194 PetscFunctionBegin; 2195 ierr = DMGlobalToLocalBegin(network->plex,g,mode,l);CHKERRQ(ierr); 2196 PetscFunctionReturn(0); 2197 } 2198 2199 PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l) 2200 { 2201 PetscErrorCode ierr; 2202 DM_Network *network = (DM_Network*)dm->data; 2203 2204 PetscFunctionBegin; 2205 ierr = DMGlobalToLocalEnd(network->plex,g,mode,l);CHKERRQ(ierr); 2206 PetscFunctionReturn(0); 2207 } 2208 2209 PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g) 2210 { 2211 PetscErrorCode ierr; 2212 DM_Network *network = (DM_Network*)dm->data; 2213 2214 PetscFunctionBegin; 2215 ierr = DMLocalToGlobalBegin(network->plex,l,mode,g);CHKERRQ(ierr); 2216 PetscFunctionReturn(0); 2217 } 2218 2219 PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g) 2220 { 2221 PetscErrorCode ierr; 2222 DM_Network *network = (DM_Network*)dm->data; 2223 2224 PetscFunctionBegin; 2225 ierr = DMLocalToGlobalEnd(network->plex,l,mode,g);CHKERRQ(ierr); 2226 PetscFunctionReturn(0); 2227 } 2228 2229 /*@ 2230 DMNetworkGetVertexLocalToGlobalOrdering - Get vertex globle index 2231 2232 Not collective 2233 2234 Input Parameters: 2235 + dm - the dm object 2236 - vloc - local vertex ordering, start from 0 2237 2238 Output Parameters: 2239 . vg - global vertex ordering, start from 0 2240 2241 Level: advanced 2242 2243 .seealso: DMNetworkSetVertexLocalToGlobalOrdering() 2244 @*/ 2245 PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg) 2246 { 2247 DM_Network *network = (DM_Network*)dm->data; 2248 PetscInt *vltog = network->vltog; 2249 2250 PetscFunctionBegin; 2251 if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first"); 2252 *vg = vltog[vloc]; 2253 PetscFunctionReturn(0); 2254 } 2255 2256 /*@ 2257 DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to globle map 2258 2259 Collective 2260 2261 Input Parameters: 2262 . dm - the dm object 2263 2264 Level: advanced 2265 2266 .seealso: DMNetworkGetGlobalVertexIndex() 2267 @*/ 2268 PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm) 2269 { 2270 PetscErrorCode ierr; 2271 DM_Network *network=(DM_Network*)dm->data; 2272 MPI_Comm comm; 2273 PetscMPIInt rank,size,*displs,*recvcounts,remoterank; 2274 PetscBool ghost; 2275 PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx; 2276 const PetscSFNode *iremote; 2277 PetscSF vsf; 2278 Vec Vleaves,Vleaves_seq; 2279 VecScatter ctx; 2280 PetscScalar *varr,val; 2281 const PetscScalar *varr_read; 2282 2283 PetscFunctionBegin; 2284 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2285 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2286 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 2287 2288 if (size == 1) { 2289 nroots = network->vEnd - network->vStart; 2290 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2291 for (i=0; i<nroots; i++) vltog[i] = i; 2292 network->vltog = vltog; 2293 PetscFunctionReturn(0); 2294 } 2295 2296 if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first"); 2297 if (network->vltog) { 2298 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2299 } 2300 2301 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 2302 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 2303 vsf = network->vertex.sf; 2304 2305 ierr = PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);CHKERRQ(ierr); 2306 ierr = PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);CHKERRQ(ierr); 2307 2308 for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;} 2309 2310 i = nroots - nleaves; /* local number of vertices, excluding ghosts */ 2311 vrange[0] = 0; 2312 ierr = MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);CHKERRQ(ierr); 2313 for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];} 2314 2315 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2316 network->vltog = vltog; 2317 2318 /* Set vltog for non-ghost vertices */ 2319 k = 0; 2320 for (i=0; i<nroots; i++) { 2321 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2322 if (ghost) continue; 2323 vltog[i] = vrange[rank] + k++; 2324 } 2325 ierr = PetscFree3(vrange,displs,recvcounts);CHKERRQ(ierr); 2326 2327 /* Set vltog for ghost vertices */ 2328 /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */ 2329 ierr = VecCreate(comm,&Vleaves);CHKERRQ(ierr); 2330 ierr = VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);CHKERRQ(ierr); 2331 ierr = VecSetFromOptions(Vleaves);CHKERRQ(ierr); 2332 ierr = VecGetArray(Vleaves,&varr);CHKERRQ(ierr); 2333 for (i=0; i<nleaves; i++) { 2334 varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */ 2335 varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */ 2336 } 2337 ierr = VecRestoreArray(Vleaves,&varr);CHKERRQ(ierr); 2338 2339 /* (b) scatter local info to remote processes via VecScatter() */ 2340 ierr = VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);CHKERRQ(ierr); 2341 ierr = VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2342 ierr = VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2343 2344 /* (c) convert local indices to global indices in parallel vector Vleaves */ 2345 ierr = VecGetSize(Vleaves_seq,&N);CHKERRQ(ierr); 2346 ierr = VecGetArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2347 for (i=0; i<N; i+=2) { 2348 remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]); 2349 if (remoterank == rank) { 2350 k = i+1; /* row number */ 2351 lidx = (PetscInt)PetscRealPart(varr_read[i+1]); 2352 val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */ 2353 ierr = VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);CHKERRQ(ierr); 2354 } 2355 } 2356 ierr = VecRestoreArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2357 ierr = VecAssemblyBegin(Vleaves);CHKERRQ(ierr); 2358 ierr = VecAssemblyEnd(Vleaves);CHKERRQ(ierr); 2359 2360 /* (d) Set vltog for ghost vertices by copying local values of Vleaves */ 2361 ierr = VecGetArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2362 k = 0; 2363 for (i=0; i<nroots; i++) { 2364 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2365 if (!ghost) continue; 2366 vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++; 2367 } 2368 ierr = VecRestoreArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2369 2370 ierr = VecDestroy(&Vleaves);CHKERRQ(ierr); 2371 ierr = VecDestroy(&Vleaves_seq);CHKERRQ(ierr); 2372 ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr); 2373 PetscFunctionReturn(0); 2374 } 2375