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