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 Parameters: 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 Parameters: 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 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: 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 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: intermediate 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: intermediate 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: intermediate 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: intermediate 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: intermediate 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: intermediate 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 /*@C 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: intermediate 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: intermediate 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: intermediate 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: intermediate 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: intermediate 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 PetscFunctionReturn(0); 1841 } 1842 1843 if (!network->userEdgeJacobian && !network->userVertexJacobian) { 1844 /* user does not provide Jacobian blocks */ 1845 ierr = DMCreateMatrix_Plex(network->plex,J);CHKERRQ(ierr); 1846 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 1847 PetscFunctionReturn(0); 1848 } 1849 1850 ierr = MatCreate(PetscObjectComm((PetscObject)dm),J);CHKERRQ(ierr); 1851 ierr = DMGetGlobalSection(network->plex,§ionGlobal);CHKERRQ(ierr); 1852 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);CHKERRQ(ierr); 1853 ierr = MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1854 1855 ierr = MatSetType(*J,MATAIJ);CHKERRQ(ierr); 1856 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 1857 1858 /* (1) Set matrix preallocation */ 1859 /*------------------------------*/ 1860 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1861 ierr = VecCreate(comm,&vd_nz);CHKERRQ(ierr); 1862 ierr = VecSetSizes(vd_nz,localSize,PETSC_DECIDE);CHKERRQ(ierr); 1863 ierr = VecSetFromOptions(vd_nz);CHKERRQ(ierr); 1864 ierr = VecSet(vd_nz,0.0);CHKERRQ(ierr); 1865 ierr = VecDuplicate(vd_nz,&vo_nz);CHKERRQ(ierr); 1866 1867 /* Set preallocation for edges */ 1868 /*-----------------------------*/ 1869 ierr = DMNetworkGetEdgeRange(dm,&eStart,&eEnd);CHKERRQ(ierr); 1870 1871 ierr = PetscMalloc1(localSize,&rows);CHKERRQ(ierr); 1872 for (e=eStart; e<eEnd; e++) { 1873 /* Get row indices */ 1874 ierr = DMNetworkGetVariableGlobalOffset(dm,e,&rstart);CHKERRQ(ierr); 1875 ierr = DMNetworkGetNumVariables(dm,e,&nrows);CHKERRQ(ierr); 1876 if (nrows) { 1877 for (j=0; j<nrows; j++) rows[j] = j + rstart; 1878 1879 /* Set preallocation for conntected vertices */ 1880 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 1881 for (v=0; v<2; v++) { 1882 ierr = DMNetworkGetNumVariables(dm,cone[v],&ncols);CHKERRQ(ierr); 1883 1884 if (network->Je) { 1885 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 1886 } else Juser = NULL; 1887 ierr = DMNetworkIsGhostVertex(dm,cone[v],&ghost);CHKERRQ(ierr); 1888 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 1889 } 1890 1891 /* Set preallocation for edge self */ 1892 cstart = rstart; 1893 if (network->Je) { 1894 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 1895 } else Juser = NULL; 1896 ierr = MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 1897 } 1898 } 1899 1900 /* Set preallocation for vertices */ 1901 /*--------------------------------*/ 1902 ierr = DMNetworkGetVertexRange(dm,&vStart,&vEnd);CHKERRQ(ierr); 1903 if (vEnd - vStart) vptr = network->Jvptr; 1904 1905 for (v=vStart; v<vEnd; v++) { 1906 /* Get row indices */ 1907 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&rstart);CHKERRQ(ierr); 1908 ierr = DMNetworkGetNumVariables(dm,v,&nrows);CHKERRQ(ierr); 1909 if (!nrows) continue; 1910 1911 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1912 if (ghost) { 1913 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 1914 } else { 1915 rows_v = rows; 1916 } 1917 1918 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 1919 1920 /* Get supporting edges and connected vertices */ 1921 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 1922 1923 for (e=0; e<nedges; e++) { 1924 /* Supporting edges */ 1925 ierr = DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);CHKERRQ(ierr); 1926 ierr = DMNetworkGetNumVariables(dm,edges[e],&ncols);CHKERRQ(ierr); 1927 1928 if (network->Jv) { 1929 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 1930 } else Juser = NULL; 1931 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);CHKERRQ(ierr); 1932 1933 /* Connected vertices */ 1934 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 1935 vc = (v == cone[0]) ? cone[1]:cone[0]; 1936 ierr = DMNetworkIsGhostVertex(dm,vc,&ghost_vc);CHKERRQ(ierr); 1937 1938 ierr = DMNetworkGetNumVariables(dm,vc,&ncols);CHKERRQ(ierr); 1939 1940 if (network->Jv) { 1941 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 1942 } else Juser = NULL; 1943 if (ghost_vc||ghost) { 1944 ghost2 = PETSC_TRUE; 1945 } else { 1946 ghost2 = PETSC_FALSE; 1947 } 1948 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);CHKERRQ(ierr); 1949 } 1950 1951 /* Set preallocation for vertex self */ 1952 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 1953 if (!ghost) { 1954 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&cstart);CHKERRQ(ierr); 1955 if (network->Jv) { 1956 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 1957 } else Juser = NULL; 1958 ierr = MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);CHKERRQ(ierr); 1959 } 1960 if (ghost) { 1961 ierr = PetscFree(rows_v);CHKERRQ(ierr); 1962 } 1963 } 1964 1965 ierr = VecAssemblyBegin(vd_nz);CHKERRQ(ierr); 1966 ierr = VecAssemblyBegin(vo_nz);CHKERRQ(ierr); 1967 1968 ierr = PetscMalloc2(localSize,&dnnz,localSize,&onnz);CHKERRQ(ierr); 1969 1970 ierr = VecAssemblyEnd(vd_nz);CHKERRQ(ierr); 1971 ierr = VecAssemblyEnd(vo_nz);CHKERRQ(ierr); 1972 1973 ierr = VecGetArray(vd_nz,&vdnz);CHKERRQ(ierr); 1974 ierr = VecGetArray(vo_nz,&vonz);CHKERRQ(ierr); 1975 for (j=0; j<localSize; j++) { 1976 dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]); 1977 onnz[j] = (PetscInt)PetscRealPart(vonz[j]); 1978 } 1979 ierr = VecRestoreArray(vd_nz,&vdnz);CHKERRQ(ierr); 1980 ierr = VecRestoreArray(vo_nz,&vonz);CHKERRQ(ierr); 1981 ierr = VecDestroy(&vd_nz);CHKERRQ(ierr); 1982 ierr = VecDestroy(&vo_nz);CHKERRQ(ierr); 1983 1984 ierr = MatSeqAIJSetPreallocation(*J,0,dnnz);CHKERRQ(ierr); 1985 ierr = MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);CHKERRQ(ierr); 1986 ierr = MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1987 1988 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 1989 1990 /* (2) Set matrix entries for edges */ 1991 /*----------------------------------*/ 1992 for (e=eStart; e<eEnd; e++) { 1993 /* Get row indices */ 1994 ierr = DMNetworkGetVariableGlobalOffset(dm,e,&rstart);CHKERRQ(ierr); 1995 ierr = DMNetworkGetNumVariables(dm,e,&nrows);CHKERRQ(ierr); 1996 if (nrows) { 1997 for (j=0; j<nrows; j++) rows[j] = j + rstart; 1998 1999 /* Set matrix entries for conntected vertices */ 2000 ierr = DMNetworkGetConnectedVertices(dm,e,&cone);CHKERRQ(ierr); 2001 for (v=0; v<2; v++) { 2002 ierr = DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);CHKERRQ(ierr); 2003 ierr = DMNetworkGetNumVariables(dm,cone[v],&ncols);CHKERRQ(ierr); 2004 2005 if (network->Je) { 2006 Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */ 2007 } else Juser = NULL; 2008 ierr = MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);CHKERRQ(ierr); 2009 } 2010 2011 /* Set matrix entries for edge self */ 2012 cstart = rstart; 2013 if (network->Je) { 2014 Juser = network->Je[3*e]; /* Jacobian(e,e) */ 2015 } else Juser = NULL; 2016 ierr = MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);CHKERRQ(ierr); 2017 } 2018 } 2019 2020 /* Set matrix entries for vertices */ 2021 /*---------------------------------*/ 2022 for (v=vStart; v<vEnd; v++) { 2023 /* Get row indices */ 2024 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&rstart);CHKERRQ(ierr); 2025 ierr = DMNetworkGetNumVariables(dm,v,&nrows);CHKERRQ(ierr); 2026 if (!nrows) continue; 2027 2028 ierr = DMNetworkIsGhostVertex(dm,v,&ghost);CHKERRQ(ierr); 2029 if (ghost) { 2030 ierr = PetscMalloc1(nrows,&rows_v);CHKERRQ(ierr); 2031 } else { 2032 rows_v = rows; 2033 } 2034 for (j=0; j<nrows; j++) rows_v[j] = j + rstart; 2035 2036 /* Get supporting edges and connected vertices */ 2037 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 2038 2039 for (e=0; e<nedges; e++) { 2040 /* Supporting edges */ 2041 ierr = DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);CHKERRQ(ierr); 2042 ierr = DMNetworkGetNumVariables(dm,edges[e],&ncols);CHKERRQ(ierr); 2043 2044 if (network->Jv) { 2045 Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */ 2046 } else Juser = NULL; 2047 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 2048 2049 /* Connected vertices */ 2050 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 2051 vc = (v == cone[0]) ? cone[1]:cone[0]; 2052 2053 ierr = DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);CHKERRQ(ierr); 2054 ierr = DMNetworkGetNumVariables(dm,vc,&ncols);CHKERRQ(ierr); 2055 2056 if (network->Jv) { 2057 Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */ 2058 } else Juser = NULL; 2059 ierr = MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);CHKERRQ(ierr); 2060 } 2061 2062 /* Set matrix entries for vertex self */ 2063 if (!ghost) { 2064 ierr = DMNetworkGetVariableGlobalOffset(dm,v,&cstart);CHKERRQ(ierr); 2065 if (network->Jv) { 2066 Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */ 2067 } else Juser = NULL; 2068 ierr = MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);CHKERRQ(ierr); 2069 } 2070 if (ghost) { 2071 ierr = PetscFree(rows_v);CHKERRQ(ierr); 2072 } 2073 } 2074 ierr = PetscFree(rows);CHKERRQ(ierr); 2075 2076 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2077 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2078 2079 ierr = MatSetDM(*J,dm);CHKERRQ(ierr); 2080 PetscFunctionReturn(0); 2081 } 2082 2083 PetscErrorCode DMDestroy_Network(DM dm) 2084 { 2085 PetscErrorCode ierr; 2086 DM_Network *network = (DM_Network*)dm->data; 2087 PetscInt j; 2088 2089 PetscFunctionBegin; 2090 if (--network->refct > 0) PetscFunctionReturn(0); 2091 if (network->Je) { 2092 ierr = PetscFree(network->Je);CHKERRQ(ierr); 2093 } 2094 if (network->Jv) { 2095 ierr = PetscFree(network->Jvptr);CHKERRQ(ierr); 2096 ierr = PetscFree(network->Jv);CHKERRQ(ierr); 2097 } 2098 2099 ierr = ISLocalToGlobalMappingDestroy(&network->vertex.mapping);CHKERRQ(ierr); 2100 ierr = PetscSectionDestroy(&network->vertex.DofSection);CHKERRQ(ierr); 2101 ierr = PetscSectionDestroy(&network->vertex.GlobalDofSection);CHKERRQ(ierr); 2102 if (network->vltog) { 2103 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2104 } 2105 if (network->vertex.sf) { 2106 ierr = PetscSFDestroy(&network->vertex.sf);CHKERRQ(ierr); 2107 } 2108 /* edge */ 2109 ierr = ISLocalToGlobalMappingDestroy(&network->edge.mapping);CHKERRQ(ierr); 2110 ierr = PetscSectionDestroy(&network->edge.DofSection);CHKERRQ(ierr); 2111 ierr = PetscSectionDestroy(&network->edge.GlobalDofSection);CHKERRQ(ierr); 2112 if (network->edge.sf) { 2113 ierr = PetscSFDestroy(&network->edge.sf);CHKERRQ(ierr); 2114 } 2115 ierr = DMDestroy(&network->plex);CHKERRQ(ierr); 2116 ierr = PetscSectionDestroy(&network->DataSection);CHKERRQ(ierr); 2117 ierr = PetscSectionDestroy(&network->DofSection);CHKERRQ(ierr); 2118 2119 for(j=0; j<network->nsubnet; j++) { 2120 ierr = PetscFree(network->subnet[j].edges);CHKERRQ(ierr); 2121 } 2122 ierr = PetscFree(network->subnetvtx);CHKERRQ(ierr); 2123 2124 ierr = PetscFree(network->subnet);CHKERRQ(ierr); 2125 ierr = PetscFree(network->componentdataarray);CHKERRQ(ierr); 2126 ierr = PetscFree2(network->header,network->cvalue);CHKERRQ(ierr); 2127 ierr = PetscFree(network);CHKERRQ(ierr); 2128 PetscFunctionReturn(0); 2129 } 2130 2131 PetscErrorCode DMView_Network(DM dm,PetscViewer viewer) 2132 { 2133 PetscErrorCode ierr; 2134 DM_Network *network = (DM_Network*)dm->data; 2135 PetscBool iascii; 2136 PetscMPIInt rank; 2137 PetscInt p,nsubnet; 2138 2139 PetscFunctionBegin; 2140 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first"); 2141 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 2142 PetscValidHeaderSpecific(dm,DM_CLASSID, 1); 2143 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 2144 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2145 if (iascii) { 2146 const PetscInt *cone,*vtx,*edges; 2147 PetscInt vfrom,vto,i,j,nv,ne; 2148 2149 nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */ 2150 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 2151 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);CHKERRQ(ierr); 2152 2153 for (i=0; i<nsubnet; i++) { 2154 ierr = DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2155 if (ne) { 2156 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);CHKERRQ(ierr); 2157 for (j=0; j<ne; j++) { 2158 p = edges[j]; 2159 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2160 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2161 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2162 ierr = DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);CHKERRQ(ierr); 2163 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2164 } 2165 } 2166 } 2167 /* Coupling subnets */ 2168 nsubnet = network->nsubnet; 2169 for (; i<nsubnet; i++) { 2170 ierr = DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 2171 if (ne) { 2172 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);CHKERRQ(ierr); 2173 for (j=0; j<ne; j++) { 2174 p = edges[j]; 2175 ierr = DMNetworkGetConnectedVertices(dm,p,&cone);CHKERRQ(ierr); 2176 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);CHKERRQ(ierr); 2177 ierr = DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);CHKERRQ(ierr); 2178 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);CHKERRQ(ierr); 2179 } 2180 } 2181 } 2182 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2183 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 2184 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name); 2185 PetscFunctionReturn(0); 2186 } 2187 2188 PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l) 2189 { 2190 PetscErrorCode ierr; 2191 DM_Network *network = (DM_Network*)dm->data; 2192 2193 PetscFunctionBegin; 2194 ierr = DMGlobalToLocalBegin(network->plex,g,mode,l);CHKERRQ(ierr); 2195 PetscFunctionReturn(0); 2196 } 2197 2198 PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l) 2199 { 2200 PetscErrorCode ierr; 2201 DM_Network *network = (DM_Network*)dm->data; 2202 2203 PetscFunctionBegin; 2204 ierr = DMGlobalToLocalEnd(network->plex,g,mode,l);CHKERRQ(ierr); 2205 PetscFunctionReturn(0); 2206 } 2207 2208 PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g) 2209 { 2210 PetscErrorCode ierr; 2211 DM_Network *network = (DM_Network*)dm->data; 2212 2213 PetscFunctionBegin; 2214 ierr = DMLocalToGlobalBegin(network->plex,l,mode,g);CHKERRQ(ierr); 2215 PetscFunctionReturn(0); 2216 } 2217 2218 PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g) 2219 { 2220 PetscErrorCode ierr; 2221 DM_Network *network = (DM_Network*)dm->data; 2222 2223 PetscFunctionBegin; 2224 ierr = DMLocalToGlobalEnd(network->plex,l,mode,g);CHKERRQ(ierr); 2225 PetscFunctionReturn(0); 2226 } 2227 2228 /*@ 2229 DMNetworkGetVertexLocalToGlobalOrdering - Get vertex globle index 2230 2231 Not collective 2232 2233 Input Parameters 2234 + dm - the dm object 2235 - vloc - local vertex ordering, start from 0 2236 2237 Output Parameters 2238 + vg - global vertex ordering, start from 0 2239 2240 Level: Advanced 2241 2242 .seealso: DMNetworkSetVertexLocalToGlobalOrdering() 2243 @*/ 2244 PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg) 2245 { 2246 DM_Network *network = (DM_Network*)dm->data; 2247 PetscInt *vltog = network->vltog; 2248 2249 PetscFunctionBegin; 2250 if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first"); 2251 *vg = vltog[vloc]; 2252 PetscFunctionReturn(0); 2253 } 2254 2255 /*@ 2256 DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to globle map 2257 2258 Collective 2259 2260 Input Parameters: 2261 + dm - the dm object 2262 2263 Level: Advanced 2264 2265 .seealso: DMNetworkGetGlobalVertexIndex() 2266 @*/ 2267 PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm) 2268 { 2269 PetscErrorCode ierr; 2270 DM_Network *network=(DM_Network*)dm->data; 2271 MPI_Comm comm; 2272 PetscMPIInt rank,size,*displs,*recvcounts,remoterank; 2273 PetscBool ghost; 2274 PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx; 2275 const PetscSFNode *iremote; 2276 PetscSF vsf; 2277 Vec Vleaves,Vleaves_seq; 2278 VecScatter ctx; 2279 PetscScalar *varr,val; 2280 const PetscScalar *varr_read; 2281 2282 PetscFunctionBegin; 2283 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2284 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2285 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 2286 2287 if (size == 1) { 2288 nroots = network->vEnd - network->vStart; 2289 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2290 for (i=0; i<nroots; i++) vltog[i] = i; 2291 network->vltog = vltog; 2292 PetscFunctionReturn(0); 2293 } 2294 2295 if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first"); 2296 if (network->vltog) { 2297 ierr = PetscFree(network->vltog);CHKERRQ(ierr); 2298 } 2299 2300 ierr = DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);CHKERRQ(ierr); 2301 ierr = PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);CHKERRQ(ierr); 2302 vsf = network->vertex.sf; 2303 2304 ierr = PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);CHKERRQ(ierr); 2305 ierr = PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);CHKERRQ(ierr); 2306 2307 for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;} 2308 2309 i = nroots - nleaves; /* local number of vertices, excluding ghosts */ 2310 vrange[0] = 0; 2311 ierr = MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);CHKERRQ(ierr); 2312 for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];} 2313 2314 ierr = PetscMalloc1(nroots, &vltog);CHKERRQ(ierr); 2315 network->vltog = vltog; 2316 2317 /* Set vltog for non-ghost vertices */ 2318 k = 0; 2319 for (i=0; i<nroots; i++) { 2320 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2321 if (ghost) continue; 2322 vltog[i] = vrange[rank] + k++; 2323 } 2324 ierr = PetscFree3(vrange,displs,recvcounts);CHKERRQ(ierr); 2325 2326 /* Set vltog for ghost vertices */ 2327 /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */ 2328 ierr = VecCreate(comm,&Vleaves);CHKERRQ(ierr); 2329 ierr = VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);CHKERRQ(ierr); 2330 ierr = VecSetFromOptions(Vleaves);CHKERRQ(ierr); 2331 ierr = VecGetArray(Vleaves,&varr);CHKERRQ(ierr); 2332 for (i=0; i<nleaves; i++) { 2333 varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */ 2334 varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */ 2335 } 2336 ierr = VecRestoreArray(Vleaves,&varr);CHKERRQ(ierr); 2337 2338 /* (b) scatter local info to remote processes via VecScatter() */ 2339 ierr = VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);CHKERRQ(ierr); 2340 ierr = VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2341 ierr = VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2342 2343 /* (c) convert local indices to global indices in parallel vector Vleaves */ 2344 ierr = VecGetSize(Vleaves_seq,&N);CHKERRQ(ierr); 2345 ierr = VecGetArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2346 for (i=0; i<N; i+=2) { 2347 remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]); 2348 if (remoterank == rank) { 2349 k = i+1; /* row number */ 2350 lidx = (PetscInt)PetscRealPart(varr_read[i+1]); 2351 val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */ 2352 ierr = VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);CHKERRQ(ierr); 2353 } 2354 } 2355 ierr = VecRestoreArrayRead(Vleaves_seq,&varr_read);CHKERRQ(ierr); 2356 ierr = VecAssemblyBegin(Vleaves);CHKERRQ(ierr); 2357 ierr = VecAssemblyEnd(Vleaves);CHKERRQ(ierr); 2358 2359 /* (d) Set vltog for ghost vertices by copying local values of Vleaves */ 2360 ierr = VecGetArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2361 k = 0; 2362 for (i=0; i<nroots; i++) { 2363 ierr = DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);CHKERRQ(ierr); 2364 if (!ghost) continue; 2365 vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++; 2366 } 2367 ierr = VecRestoreArrayRead(Vleaves,&varr_read);CHKERRQ(ierr); 2368 2369 ierr = VecDestroy(&Vleaves);CHKERRQ(ierr); 2370 ierr = VecDestroy(&Vleaves_seq);CHKERRQ(ierr); 2371 ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr); 2372 PetscFunctionReturn(0); 2373 } 2374