Actual source code: network.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/dmnetworkimpl.h>
2: #include <petscdmplex.h>
3: #include <petscsf.h>
5: /*@
6: DMNetworkGetPlex - Gets the Plex DM associated with this network DM
8: Not collective
10: Input Parameters:
11: + netdm - the dm object
12: - plexmdm - the plex dm object
14: Level: Advanced
16: .seealso: DMNetworkCreate()
17: @*/
18: PetscErrorCode DMNetworkGetPlex(DM netdm, DM *plexdm)
19: {
20: DM_Network *network = (DM_Network*) netdm->data;
23: *plexdm = network->plex;
24: return(0);
25: }
27: /*@
28: DMNetworkSetSizes - Sets the local and global vertices and edges.
30: Collective on DM
32: Input Parameters:
33: + dm - the dm object
34: . nV - number of local vertices
35: . nE - number of local edges
36: . NV - number of global vertices (or PETSC_DETERMINE)
37: - NE - number of global edges (or PETSC_DETERMINE)
39: Notes
40: If one processor calls this with NV (NE) of PETSC_DECIDE then all processors must, otherwise the prgram will hang.
42: You cannot change the sizes once they have been set
44: Level: intermediate
46: .seealso: DMNetworkCreate()
47: @*/
48: PetscErrorCode DMNetworkSetSizes(DM dm, PetscInt nV, PetscInt nE, PetscInt NV, PetscInt NE)
49: {
51: DM_Network *network = (DM_Network*) dm->data;
52: PetscInt a[2],b[2];
58: if (NV > 0 && nV > NV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local vertex size %D cannot be larger than global vertex size %D",nV,NV);
59: if (NE > 0 && nE > NE) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local edge size %D cannot be larger than global edge size %D",nE,NE);
60: if ((network->nVertices >= 0 || network->NVertices >= 0) && (network->nVertices != nV || network->NVertices != NV)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset vertex sizes to %D local %D global after previously setting them to %D local %D global",nV,NV,network->nVertices,network->NVertices);
61: if ((network->nEdges >= 0 || network->NEdges >= 0) && (network->nEdges != nE || network->NEdges != NE)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset edge sizes to %D local %D global after previously setting them to %D local %D global",nE,NE,network->nEdges,network->NEdges);
62: if (NE < 0 || NV < 0) {
63: a[0] = nV; a[1] = nE;
64: MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
65: NV = b[0]; NE = b[1];
66: }
67: network->nVertices = nV;
68: network->NVertices = NV;
69: network->nEdges = nE;
70: network->NEdges = NE;
71: return(0);
72: }
74: /*@
75: DMNetworkSetEdgeList - Sets the list of local edges (vertex connectivity) for the network
77: Logically collective on DM
79: Input Parameters:
80: . edges - list of edges
82: Notes:
83: There is no copy involved in this operation, only the pointer is referenced. The edgelist should
84: not be destroyed before the call to DMNetworkLayoutSetUp
86: Level: intermediate
88: .seealso: DMNetworkCreate, DMNetworkSetSizes
89: @*/
90: PetscErrorCode DMNetworkSetEdgeList(DM dm, int edgelist[])
91: {
92: DM_Network *network = (DM_Network*) dm->data;
95: network->edges = edgelist;
96: return(0);
97: }
99: /*@
100: DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network
102: Collective on DM
104: Input Parameters
105: . DM - the dmnetwork object
107: Notes:
108: This routine should be called after the network sizes and edgelists have been provided. It creates
109: the bare layout of the network and sets up the network to begin insertion of components.
111: All the components should be registered before calling this routine.
113: Level: intermediate
115: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
116: @*/
117: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
118: {
120: DM_Network *network = (DM_Network*) dm->data;
121: PetscInt dim = 1; /* One dimensional network */
122: PetscInt numCorners=2;
123: PetscInt spacedim=2;
124: double *vertexcoords=NULL;
125: PetscInt i;
126: PetscInt ndata;
129: if (network->nVertices) {
130: PetscCalloc1(numCorners*network->nVertices,&vertexcoords);
131: }
132: DMPlexCreateFromCellList(PetscObjectComm((PetscObject)dm),dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,network->edges,spacedim,vertexcoords,&network->plex);
133: if (network->nVertices) {
134: PetscFree(vertexcoords);
135: }
136: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
137: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
138: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
140: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DataSection);
141: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DofSection);
142: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
143: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
147: network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
148: PetscCalloc1(network->pEnd-network->pStart,&network->header);
149: for (i = network->pStart; i < network->pEnd; i++) {
150: if (i < network->vStart) {
151: network->header[i].index = i;
152: } else {
153: network->header[i].index = i - network->vStart;
154: }
155: network->header[i].ndata = 0;
156: ndata = network->header[i].ndata;
157: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
158: network->header[i].offset[ndata] = 0;
159: }
160: PetscMalloc1(network->pEnd-network->pStart,&network->cvalue);
161: return(0);
162: }
164: /*@C
165: DMNetworkRegisterComponent - Registers the network component
167: Logically collective on DM
169: Input Parameters
170: + dm - the network object
171: . name - the component name
172: - size - the storage size in bytes for this component data
174: Output Parameters
175: . key - an integer key that defines the component
177: Notes
178: This routine should be called by all processors before calling DMNetworkLayoutSetup().
180: Level: intermediate
182: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
183: @*/
184: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,PetscInt size,PetscInt *key)
185: {
186: PetscErrorCode ierr;
187: DM_Network *network = (DM_Network*) dm->data;
188: DMNetworkComponent *component=&network->component[network->ncomponent];
189: PetscBool flg=PETSC_FALSE;
190: PetscInt i;
194: for (i=0; i < network->ncomponent; i++) {
195: PetscStrcmp(component->name,name,&flg);
196: if (flg) {
197: *key = i;
198: return(0);
199: }
200: }
202: PetscStrcpy(component->name,name);
203: component->size = size/sizeof(DMNetworkComponentGenericDataType);
204: *key = network->ncomponent;
205: network->ncomponent++;
206: return(0);
207: }
209: /*@
210: DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices.
212: Not Collective
214: Input Parameters:
215: + dm - The DMNetwork object
217: Output Paramters:
218: + vStart - The first vertex point
219: - vEnd - One beyond the last vertex point
221: Level: intermediate
223: .seealso: DMNetworkGetEdgeRange
224: @*/
225: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
226: {
227: DM_Network *network = (DM_Network*)dm->data;
230: if (vStart) *vStart = network->vStart;
231: if (vEnd) *vEnd = network->vEnd;
232: return(0);
233: }
235: /*@
236: DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges.
238: Not Collective
240: Input Parameters:
241: + dm - The DMNetwork object
243: Output Paramters:
244: + eStart - The first edge point
245: - eEnd - One beyond the last edge point
247: Level: intermediate
249: .seealso: DMNetworkGetVertexRange
250: @*/
251: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
252: {
253: DM_Network *network = (DM_Network*)dm->data;
256: if (eStart) *eStart = network->eStart;
257: if (eEnd) *eEnd = network->eEnd;
258: return(0);
259: }
261: /*@
262: DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.
264: Not Collective
266: Input Parameters:
267: + dm - DMNetwork object
268: - p - edge point
270: Output Paramters:
271: . index - user global numbering for the edge
273: Level: intermediate
275: .seealso: DMNetworkGetGlobalVertexIndex
276: @*/
277: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
278: {
279: PetscErrorCode ierr;
280: DM_Network *network = (DM_Network*)dm->data;
281: PetscInt offsetp;
282: DMNetworkComponentHeader header;
285: PetscSectionGetOffset(network->DataSection,p,&offsetp);
286: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
287: *index = header->index;
288: return(0);
289: }
291: /*@
292: DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.
294: Not Collective
296: Input Parameters:
297: + dm - DMNetwork object
298: - p - vertex point
300: Output Paramters:
301: . index - user global numbering for the vertex
303: Level: intermediate
305: .seealso: DMNetworkGetGlobalEdgeIndex
306: @*/
307: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
308: {
309: PetscErrorCode ierr;
310: DM_Network *network = (DM_Network*)dm->data;
311: PetscInt offsetp;
312: DMNetworkComponentHeader header;
315: PetscSectionGetOffset(network->DataSection,p,&offsetp);
316: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
317: *index = header->index;
318: return(0);
319: }
321: /*@
322: DMNetworkAddComponent - Adds a network component at the given point (vertex/edge)
324: Not Collective
326: Input Parameters:
327: + dm - The DMNetwork object
328: . p - vertex/edge point
329: . componentkey - component key returned while registering the component
330: - compvalue - pointer to the data structure for the component
332: Level: intermediate
334: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
335: @*/
336: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
337: {
338: DM_Network *network = (DM_Network*)dm->data;
339: DMNetworkComponent *component = &network->component[componentkey];
340: DMNetworkComponentHeader header = &network->header[p];
341: DMNetworkComponentValue cvalue = &network->cvalue[p];
342: PetscErrorCode ierr;
345: 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);
347: header->size[header->ndata] = component->size;
348: PetscSectionAddDof(network->DataSection,p,component->size);
349: header->key[header->ndata] = componentkey;
350: if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];
352: cvalue->data[header->ndata] = (void*)compvalue;
353: header->ndata++;
354: return(0);
355: }
357: /*@
358: DMNetworkGetNumComponents - Get the number of components at a vertex/edge
360: Not Collective
362: Input Parameters:
363: + dm - The DMNetwork object
364: . p - vertex/edge point
366: Output Parameters:
367: . numcomponents - Number of components at the vertex/edge
369: Level: intermediate
371: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
372: @*/
373: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
374: {
376: PetscInt offset;
377: DM_Network *network = (DM_Network*)dm->data;
380: PetscSectionGetOffset(network->DataSection,p,&offset);
381: *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
382: return(0);
383: }
385: /*@
386: DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
387: component value from the component data array
389: Not Collective
391: Input Parameters:
392: + dm - The DMNetwork object
393: . p - vertex/edge point
394: - compnum - component number
396: Output Parameters:
397: + compkey - the key obtained when registering the component
398: - offset - offset into the component data array associated with the vertex/edge point
400: Notes:
401: Typical usage:
403: DMNetworkGetComponentDataArray(dm, &arr);
404: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
405: Loop over vertices or edges
406: DMNetworkGetNumComponents(dm,v,&numcomps);
407: Loop over numcomps
408: DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
409: compdata = (UserCompDataType)(arr+offset);
411: Level: intermediate
413: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
414: @*/
415: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
416: {
417: PetscErrorCode ierr;
418: PetscInt offsetp;
419: DMNetworkComponentHeader header;
420: DM_Network *network = (DM_Network*)dm->data;
423: PetscSectionGetOffset(network->DataSection,p,&offsetp);
424: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
425: if (compkey) *compkey = header->key[compnum];
426: if (offset) *offset = offsetp+network->dataheadersize+header->offset[compnum];
427: return(0);
428: }
430: /*@
431: DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector.
433: Not Collective
435: Input Parameters:
436: + dm - The DMNetwork object
437: - p - the edge/vertex point
439: Output Parameters:
440: . offset - the offset
442: Level: intermediate
444: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
445: @*/
446: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
447: {
449: DM_Network *network = (DM_Network*)dm->data;
452: PetscSectionGetOffset(network->plex->defaultSection,p,offset);
453: return(0);
454: }
456: /*@
457: DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector.
459: Not Collective
461: Input Parameters:
462: + dm - The DMNetwork object
463: - p - the edge/vertex point
465: Output Parameters:
466: . offsetg - the offset
468: Level: intermediate
470: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
471: @*/
472: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
473: {
475: DM_Network *network = (DM_Network*)dm->data;
478: PetscSectionGetOffset(network->plex->defaultGlobalSection,p,offsetg);
479: if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
480: return(0);
481: }
483: /*@
484: DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector.
486: Not Collective
488: Input Parameters:
489: + dm - The DMNetwork object
490: - p - the edge point
492: Output Parameters:
493: . offset - the offset
495: Level: intermediate
497: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
498: @*/
499: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
500: {
502: DM_Network *network = (DM_Network*)dm->data;
506: PetscSectionGetOffset(network->edge.DofSection,p,offset);
507: return(0);
508: }
510: /*@
511: DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector.
513: Not Collective
515: Input Parameters:
516: + dm - The DMNetwork object
517: - p - the vertex point
519: Output Parameters:
520: . offset - the offset
522: Level: intermediate
524: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
525: @*/
526: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
527: {
529: DM_Network *network = (DM_Network*)dm->data;
533: p -= network->vStart;
535: PetscSectionGetOffset(network->vertex.DofSection,p,offset);
536: return(0);
537: }
538: /*@
539: DMNetworkAddNumVariables - Add number of variables associated with a given point.
541: Not Collective
543: Input Parameters:
544: + dm - The DMNetworkObject
545: . p - the vertex/edge point
546: - nvar - number of additional variables
548: Level: intermediate
550: .seealso: DMNetworkSetNumVariables
551: @*/
552: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
553: {
555: DM_Network *network = (DM_Network*)dm->data;
558: PetscSectionAddDof(network->DofSection,p,nvar);
559: return(0);
560: }
562: /*@
563: DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point.
565: Not Collective
567: Input Parameters:
568: + dm - The DMNetworkObject
569: - p - the vertex/edge point
571: Output Parameters:
572: . nvar - number of variables
574: Level: intermediate
576: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
577: @*/
578: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
579: {
581: DM_Network *network = (DM_Network*)dm->data;
584: PetscSectionGetDof(network->DofSection,p,nvar);
585: return(0);
586: }
588: /*@
589: DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point.
591: Not Collective
593: Input Parameters:
594: + dm - The DMNetworkObject
595: . p - the vertex/edge point
596: - nvar - number of variables
598: Level: intermediate
600: .seealso: DMNetworkAddNumVariables
601: @*/
602: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
603: {
605: DM_Network *network = (DM_Network*)dm->data;
608: PetscSectionSetDof(network->DofSection,p,nvar);
609: return(0);
610: }
612: /* Sets up the array that holds the data for all components and its associated section. This
613: function is called during DMSetUp() */
614: PetscErrorCode DMNetworkComponentSetUp(DM dm)
615: {
616: PetscErrorCode ierr;
617: DM_Network *network = (DM_Network*)dm->data;
618: PetscInt arr_size;
619: PetscInt p,offset,offsetp;
620: DMNetworkComponentHeader header;
621: DMNetworkComponentValue cvalue;
622: DMNetworkComponentGenericDataType *componentdataarray;
623: PetscInt ncomp, i;
626: PetscSectionSetUp(network->DataSection);
627: PetscSectionGetStorageSize(network->DataSection,&arr_size);
628: PetscMalloc1(arr_size,&network->componentdataarray);
629: componentdataarray = network->componentdataarray;
630: for (p = network->pStart; p < network->pEnd; p++) {
631: PetscSectionGetOffset(network->DataSection,p,&offsetp);
632: /* Copy header */
633: header = &network->header[p];
634: PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
635: /* Copy data */
636: cvalue = &network->cvalue[p];
637: ncomp = header->ndata;
638: for (i = 0; i < ncomp; i++) {
639: offset = offsetp + network->dataheadersize + header->offset[i];
640: PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
641: }
642: }
643: return(0);
644: }
646: /* Sets up the section for dofs. This routine is called during DMSetUp() */
647: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
648: {
650: DM_Network *network = (DM_Network*)dm->data;
653: PetscSectionSetUp(network->DofSection);
654: return(0);
655: }
657: /*@C
658: DMNetworkGetComponentDataArray - Returns the component data array
660: Not Collective
662: Input Parameters:
663: . dm - The DMNetwork Object
665: Output Parameters:
666: . componentdataarray - array that holds data for all components
668: Level: intermediate
670: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
671: @*/
672: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
673: {
674: DM_Network *network = (DM_Network*)dm->data;
677: *componentdataarray = network->componentdataarray;
678: return(0);
679: }
681: /* Get a subsection from a range of points */
682: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
683: {
685: PetscInt i, nvar;
688: PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
689: PetscSectionSetChart(*subsection, 0, pend - pstart);
690: for (i = pstart; i < pend; i++) {
691: PetscSectionGetDof(master,i,&nvar);
692: PetscSectionSetDof(*subsection, i - pstart, nvar);
693: }
695: PetscSectionSetUp(*subsection);
696: return(0);
697: }
699: /* Create a submap of points with a GlobalToLocal structure */
700: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
701: {
703: PetscInt i, *subpoints;
706: /* Create index sets to map from "points" to "subpoints" */
707: PetscMalloc1(pend - pstart, &subpoints);
708: for (i = pstart; i < pend; i++) {
709: subpoints[i - pstart] = i;
710: }
711: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
712: PetscFree(subpoints);
713: return(0);
714: }
716: /*@
717: DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute.
719: Collective
721: Input Parameters:
722: . dm - The DMNetworkObject
724: Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are:
726: points = [0 1 2 3 4 5 6]
728: where edges = [0, 3] and vertices = [4, 6]. The new orderings will be specific to the subset (i.e vertices = [0, 2]).
730: With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset.
732: Level: intermediate
734: @*/
735: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
736: {
738: MPI_Comm comm;
739: PetscMPIInt rank, size;
740: DM_Network *network = (DM_Network*)dm->data;
743: PetscObjectGetComm((PetscObject)dm,&comm);
744: MPI_Comm_rank(comm, &rank);
745: MPI_Comm_size(comm, &size);
747: /* Create maps for vertices and edges */
748: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
749: DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);
751: /* Create local sub-sections */
752: DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);
753: DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);
755: if (size > 1) {
756: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
757: PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
758: PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
759: PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
760: } else {
761: /* create structures for vertex */
762: PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
763: /* create structures for edge */
764: PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
765: }
768: /* Add viewers */
769: PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
770: PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
771: PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
772: PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
774: return(0);
775: }
777: /*@
778: DMNetworkDistribute - Distributes the network and moves associated component data.
780: Collective
782: Input Parameter:
783: + DM - the DMNetwork object
784: - overlap - The overlap of partitions, 0 is the default
786: Notes:
787: Distributes the network with <overlap>-overlapping partitioning of the edges.
789: Level: intermediate
791: .seealso: DMNetworkCreate
792: @*/
793: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
794: {
795: MPI_Comm comm;
797: PetscMPIInt size;
798: DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data);
799: DM_Network *newDMnetwork;
800: PetscSF pointsf;
801: DM newDM;
802: PetscPartitioner part;
806: PetscObjectGetComm((PetscObject)*dm,&comm);
807: MPI_Comm_size(comm, &size);
808: if (size == 1) return(0);
810: DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
811: newDMnetwork = (DM_Network*)newDM->data;
812: newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
814: /* Enable runtime options for petscpartitioner */
815: DMPlexGetPartitioner(oldDMnetwork->plex,&part);
816: PetscPartitionerSetFromOptions(part);
818: /* Distribute plex dm and dof section */
819: DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
821: /* Distribute dof section */
822: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
823: PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
824: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);
826: /* Distribute data and associated section */
827: DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
829: PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
830: DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
831: DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
832: newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
833: newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
834: newDMnetwork->NVertices = oldDMnetwork->NVertices;
835: newDMnetwork->NEdges = oldDMnetwork->NEdges;
837: /* Set Dof section as the default section for dm */
838: DMSetDefaultSection(newDMnetwork->plex,newDMnetwork->DofSection);
839: DMGetDefaultGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);
841: /* Destroy point SF */
842: PetscSFDestroy(&pointsf);
844: DMDestroy(dm);
845: *dm = newDM;
846: return(0);
847: }
849: /*@C
850: PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering.
852: Input Parameters:
853: + masterSF - the original SF structure
854: - map - a ISLocalToGlobal mapping that contains the subset of points
856: Output Parameters:
857: . subSF - a subset of the masterSF for the desired subset.
858: */
860: PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) {
862: PetscErrorCode ierr;
863: PetscInt nroots, nleaves, *ilocal_sub;
864: PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0;
865: PetscInt *local_points, *remote_points;
866: PetscSFNode *iremote_sub;
867: const PetscInt *ilocal;
868: const PetscSFNode *iremote;
871: PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);
873: /* Look for leaves that pertain to the subset of points. Get the local ordering */
874: PetscMalloc1(nleaves,&ilocal_map);
875: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
876: for (i = 0; i < nleaves; i++) {
877: if (ilocal_map[i] != -1) nleaves_sub += 1;
878: }
879: /* Re-number ilocal with subset numbering. Need information from roots */
880: PetscMalloc2(nroots,&local_points,nroots,&remote_points);
881: for (i = 0; i < nroots; i++) local_points[i] = i;
882: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
883: PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
884: PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
885: /* Fill up graph using local (that is, local to the subset) numbering. */
886: PetscMalloc1(nleaves_sub,&ilocal_sub);
887: PetscMalloc1(nleaves_sub,&iremote_sub);
888: nleaves_sub = 0;
889: for (i = 0; i < nleaves; i++) {
890: if (ilocal_map[i] != -1) {
891: ilocal_sub[nleaves_sub] = ilocal_map[i];
892: iremote_sub[nleaves_sub].rank = iremote[i].rank;
893: iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
894: nleaves_sub += 1;
895: }
896: }
897: PetscFree2(local_points,remote_points);
898: ISLocalToGlobalMappingGetSize(map,&nroots_sub);
900: /* Create new subSF */
901: PetscSFCreate(PETSC_COMM_WORLD,subSF);
902: PetscSFSetFromOptions(*subSF);
903: PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
904: PetscFree(ilocal_map);
905: PetscFree(iremote_sub);
906: return(0);
907: }
909: /*@C
910: DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point
912: Not Collective
914: Input Parameters:
915: + dm - The DMNetwork object
916: - p - the vertex point
918: Output Paramters:
919: + nedges - number of edges connected to this vertex point
920: - edges - List of edge points
922: Level: intermediate
924: Fortran Notes:
925: Since it returns an array, this routine is only available in Fortran 90, and you must
926: include petsc.h90 in your code.
928: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
929: @*/
930: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
931: {
933: DM_Network *network = (DM_Network*)dm->data;
936: DMPlexGetSupportSize(network->plex,vertex,nedges);
937: DMPlexGetSupport(network->plex,vertex,edges);
938: return(0);
939: }
941: /*@C
942: DMNetworkGetConnectedVertices - Return the connected vertices for this edge point
944: Not Collective
946: Input Parameters:
947: + dm - The DMNetwork object
948: - p - the edge point
950: Output Paramters:
951: . vertices - vertices connected to this edge
953: Level: intermediate
955: Fortran Notes:
956: Since it returns an array, this routine is only available in Fortran 90, and you must
957: include petsc.h90 in your code.
959: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
960: @*/
961: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
962: {
964: DM_Network *network = (DM_Network*)dm->data;
967: DMPlexGetCone(network->plex,edge,vertices);
968: return(0);
969: }
971: /*@
972: DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex
974: Not Collective
976: Input Parameters:
977: + dm - The DMNetwork object
978: . p - the vertex point
980: Output Parameter:
981: . isghost - TRUE if the vertex is a ghost point
983: Level: intermediate
985: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
986: @*/
987: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
988: {
990: DM_Network *network = (DM_Network*)dm->data;
991: PetscInt offsetg;
992: PetscSection sectiong;
995: *isghost = PETSC_FALSE;
996: DMGetDefaultGlobalSection(network->plex,§iong);
997: PetscSectionGetOffset(sectiong,p,&offsetg);
998: if (offsetg < 0) *isghost = PETSC_TRUE;
999: return(0);
1000: }
1002: PetscErrorCode DMSetUp_Network(DM dm)
1003: {
1005: DM_Network *network=(DM_Network*)dm->data;
1008: DMNetworkComponentSetUp(dm);
1009: DMNetworkVariablesSetUp(dm);
1011: DMSetDefaultSection(network->plex,network->DofSection);
1012: DMGetDefaultGlobalSection(network->plex,&network->GlobalDofSection);
1013: return(0);
1014: }
1016: /*@
1017: DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
1018: -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?
1020: Collective
1022: Input Parameters:
1023: + dm - The DMNetwork object
1024: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1025: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices
1027: Level: intermediate
1029: @*/
1030: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1031: {
1032: DM_Network *network=(DM_Network*)dm->data;
1035: network->userEdgeJacobian = eflg;
1036: network->userVertexJacobian = vflg;
1037: return(0);
1038: }
1040: /*@
1041: DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network
1043: Not Collective
1045: Input Parameters:
1046: + dm - The DMNetwork object
1047: . p - the edge point
1048: - J - array (size = 3) of Jacobian submatrices for this edge point:
1049: J[0]: this edge
1050: J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()
1052: Level: intermediate
1054: .seealso: DMNetworkVertexSetMatrix
1055: @*/
1056: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1057: {
1059: DM_Network *network=(DM_Network*)dm->data;
1062: if (!network->userEdgeJacobian) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
1063: if (!network->Je) {
1064: PetscCalloc1(3*network->nEdges,&network->Je);
1065: }
1066: network->Je[3*p] = J[0];
1067: network->Je[3*p+1] = J[1];
1068: network->Je[3*p+2] = J[2];
1069: return(0);
1070: }
1072: /*@
1073: DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network
1075: Not Collective
1077: Input Parameters:
1078: + dm - The DMNetwork object
1079: . p - the vertex point
1080: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1081: J[0]: this vertex
1082: J[1+2*i]: i-th supporting edge
1083: J[1+2*i+1]: i-th connected vertex
1085: Level: intermediate
1087: .seealso: DMNetworkEdgeSetMatrix
1088: @*/
1089: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1090: {
1092: DM_Network *network=(DM_Network*)dm->data;
1093: PetscInt i,*vptr,nedges,vStart,vEnd;
1094: const PetscInt *edges;
1097: if (!network->userVertexJacobian) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");
1099: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1101: if (!network->Jv) {
1102: PetscInt nVertices = network->nVertices,nedges_total;
1103: /* count nvertex_total */
1104: nedges_total = 0;
1105: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1106: PetscMalloc1(nVertices+1,&vptr);
1108: vptr[0] = 0;
1109: for (i=0; i<nVertices; i++) {
1110: DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1111: nedges_total += nedges;
1112: vptr[i+1] = vptr[i] + 2*nedges + 1;
1113: }
1115: PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1116: network->Jvptr = vptr;
1117: }
1119: vptr = network->Jvptr;
1120: network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */
1122: /* Set Jacobian for each supporting edge and connected vertex */
1123: DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1124: for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1125: return(0);
1126: }
1128: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1129: {
1131: PetscInt j;
1132: PetscScalar val=(PetscScalar)ncols;
1135: if (!ghost) {
1136: for (j=0; j<nrows; j++) {
1137: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1138: }
1139: } else {
1140: for (j=0; j<nrows; j++) {
1141: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1142: }
1143: }
1144: return(0);
1145: }
1147: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1148: {
1150: PetscInt j,ncols_u;
1151: PetscScalar val;
1154: if (!ghost) {
1155: for (j=0; j<nrows; j++) {
1156: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1157: val = (PetscScalar)ncols_u;
1158: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1159: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1160: }
1161: } else {
1162: for (j=0; j<nrows; j++) {
1163: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1164: val = (PetscScalar)ncols_u;
1165: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1166: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1167: }
1168: }
1169: return(0);
1170: }
1172: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1173: {
1177: if (Ju) {
1178: MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1179: } else {
1180: MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1181: }
1182: return(0);
1183: }
1185: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1186: {
1188: PetscInt j,*cols;
1189: PetscScalar *zeros;
1192: PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1193: for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1194: MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1195: PetscFree2(cols,zeros);
1196: return(0);
1197: }
1199: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1200: {
1202: PetscInt j,M,N,row,col,ncols_u;
1203: const PetscInt *cols;
1204: PetscScalar zero=0.0;
1207: MatGetSize(Ju,&M,&N);
1208: if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);
1210: for (row=0; row<nrows; row++) {
1211: MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1212: for (j=0; j<ncols_u; j++) {
1213: col = cols[j] + cstart;
1214: MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1215: }
1216: MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1217: }
1218: return(0);
1219: }
1221: PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1222: {
1226: if (Ju) {
1227: MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1228: } else {
1229: MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1230: }
1231: return(0);
1232: }
1234: /* 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.
1235: */
1236: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1237: {
1239: PetscInt i, size, dof;
1240: PetscInt *glob2loc;
1243: PetscSectionGetStorageSize(localsec,&size);
1244: PetscMalloc1(size,&glob2loc);
1246: for (i = 0; i < size; i++) {
1247: PetscSectionGetOffset(globalsec,i,&dof);
1248: dof = (dof >= 0) ? dof : -(dof + 1);
1249: glob2loc[i] = dof;
1250: }
1252: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1253: #if 0
1254: PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1255: #endif
1256: return(0);
1257: }
1259: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1260: {
1262: PetscMPIInt rank, size;
1263: DM_Network *network = (DM_Network*) dm->data;
1264: PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1265: PetscInt cstart,ncols,j,e,v;
1266: PetscBool ghost,ghost_vc,ghost2,isNest;
1267: Mat Juser;
1268: PetscSection sectionGlobal;
1269: PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1270: const PetscInt *edges,*cone;
1271: MPI_Comm comm;
1272: MatType mtype;
1273: Vec vd_nz,vo_nz;
1274: PetscInt *dnnz,*onnz;
1275: PetscScalar *vdnz,*vonz;
1278: mtype = dm->mattype;
1279: PetscStrcmp(mtype, MATNEST, &isNest);
1281: if (isNest) {
1282: /* DMCreateMatrix_Network_Nest(); */
1283: PetscInt eDof, vDof;
1284: Mat j11, j12, j21, j22, bA[2][2];
1285: ISLocalToGlobalMapping eISMap, vISMap;
1287: PetscObjectGetComm((PetscObject)dm,&comm);
1288: MPI_Comm_rank(comm,&rank);
1289: MPI_Comm_size(comm,&size);
1291: PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);
1292: PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);
1294: MatCreate(PETSC_COMM_WORLD, &j11);
1295: MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1296: MatSetType(j11, MATMPIAIJ);
1298: MatCreate(PETSC_COMM_WORLD, &j12);
1299: MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1300: MatSetType(j12, MATMPIAIJ);
1302: MatCreate(PETSC_COMM_WORLD, &j21);
1303: MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1304: MatSetType(j21, MATMPIAIJ);
1306: MatCreate(PETSC_COMM_WORLD, &j22);
1307: MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1308: MatSetType(j22, MATMPIAIJ);
1310: bA[0][0] = j11;
1311: bA[0][1] = j12;
1312: bA[1][0] = j21;
1313: bA[1][1] = j22;
1315: CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);
1316: CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);
1318: MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1319: MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1320: MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1321: MatSetLocalToGlobalMapping(j22,vISMap,vISMap);
1323: MatSetUp(j11);
1324: MatSetUp(j12);
1325: MatSetUp(j21);
1326: MatSetUp(j22);
1328: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],J);
1329: MatSetUp(*J);
1330: MatNestSetVecType(*J,VECNEST);
1331: MatDestroy(&j11);
1332: MatDestroy(&j12);
1333: MatDestroy(&j21);
1334: MatDestroy(&j22);
1336: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1337: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1339: /* Free structures */
1340: ISLocalToGlobalMappingDestroy(&eISMap);
1341: ISLocalToGlobalMappingDestroy(&vISMap);
1343: return(0);
1344: } else if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1345: /* user does not provide Jacobian blocks */
1346: DMCreateMatrix(network->plex,J);
1347: MatSetDM(*J,dm);
1348: return(0);
1349: }
1351: MatCreate(PetscObjectComm((PetscObject)dm),J);
1352: DMGetDefaultGlobalSection(network->plex,§ionGlobal);
1353: PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1354: MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
1356: MatSetType(*J,MATAIJ);
1357: MatSetFromOptions(*J);
1359: /* (1) Set matrix preallocation */
1360: /*------------------------------*/
1361: PetscObjectGetComm((PetscObject)dm,&comm);
1362: VecCreate(comm,&vd_nz);
1363: VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1364: VecSetFromOptions(vd_nz);
1365: VecSet(vd_nz,0.0);
1366: VecDuplicate(vd_nz,&vo_nz);
1368: /* Set preallocation for edges */
1369: /*-----------------------------*/
1370: DMNetworkGetEdgeRange(dm,&eStart,&eEnd);
1372: PetscMalloc1(localSize,&rows);
1373: for (e=eStart; e<eEnd; e++) {
1374: /* Get row indices */
1375: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1376: DMNetworkGetNumVariables(dm,e,&nrows);
1377: if (nrows) {
1378: if (!network->Je) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Je");
1379: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1381: /* Set preallocation for conntected vertices */
1382: DMNetworkGetConnectedVertices(dm,e,&cone);
1383: for (v=0; v<2; v++) {
1384: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1386: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1387: DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1388: MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1389: }
1391: /* Set preallocation for edge self */
1392: cstart = rstart;
1393: Juser = network->Je[3*e]; /* Jacobian(e,e) */
1394: MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1395: }
1396: }
1398: /* Set preallocation for vertices */
1399: /*--------------------------------*/
1400: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1401: if (vEnd - vStart) {
1402: if (!network->Jv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Jv");
1403: vptr = network->Jvptr;
1404: if (!vptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide vptr");
1405: }
1407: for (v=vStart; v<vEnd; v++) {
1408: /* Get row indices */
1409: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1410: DMNetworkGetNumVariables(dm,v,&nrows);
1411: if (!nrows) continue;
1413: DMNetworkIsGhostVertex(dm,v,&ghost);
1414: if (ghost) {
1415: PetscMalloc1(nrows,&rows_v);
1416: } else {
1417: rows_v = rows;
1418: }
1420: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
1422: /* Get supporting edges and connected vertices */
1423: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1425: for (e=0; e<nedges; e++) {
1426: /* Supporting edges */
1427: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1428: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1430: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1431: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);
1433: /* Connected vertices */
1434: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1435: vc = (v == cone[0]) ? cone[1]:cone[0];
1436: DMNetworkIsGhostVertex(dm,vc,&ghost_vc);
1438: DMNetworkGetNumVariables(dm,vc,&ncols);
1440: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1441: if (ghost_vc||ghost) {
1442: ghost2 = PETSC_TRUE;
1443: } else {
1444: ghost2 = PETSC_FALSE;
1445: }
1446: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1447: }
1449: /* Set preallocation for vertex self */
1450: DMNetworkIsGhostVertex(dm,v,&ghost);
1451: if (!ghost) {
1452: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1453: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1454: MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1455: }
1456: if (ghost) {
1457: PetscFree(rows_v);
1458: }
1459: }
1461: VecAssemblyBegin(vd_nz);
1462: VecAssemblyBegin(vo_nz);
1464: PetscMalloc2(localSize,&dnnz,localSize,&onnz);
1466: VecAssemblyEnd(vd_nz);
1467: VecAssemblyEnd(vo_nz);
1469: VecGetArray(vd_nz,&vdnz);
1470: VecGetArray(vo_nz,&vonz);
1471: for (j=0; j<localSize; j++) {
1472: dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1473: onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1474: }
1475: VecRestoreArray(vd_nz,&vdnz);
1476: VecRestoreArray(vo_nz,&vonz);
1477: VecDestroy(&vd_nz);
1478: VecDestroy(&vo_nz);
1480: MatSeqAIJSetPreallocation(*J,0,dnnz);
1481: MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1482: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1484: PetscFree2(dnnz,onnz);
1486: /* (2) Set matrix entries for edges */
1487: /*----------------------------------*/
1488: for (e=eStart; e<eEnd; e++) {
1489: /* Get row indices */
1490: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1491: DMNetworkGetNumVariables(dm,e,&nrows);
1492: if (nrows) {
1493: if (!network->Je) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Je");
1495: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1497: /* Set matrix entries for conntected vertices */
1498: DMNetworkGetConnectedVertices(dm,e,&cone);
1499: for (v=0; v<2; v++) {
1500: DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
1501: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1503: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1504: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
1505: }
1507: /* Set matrix entries for edge self */
1508: cstart = rstart;
1509: Juser = network->Je[3*e]; /* Jacobian(e,e) */
1510: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1511: }
1512: }
1514: /* Set matrix entries for vertices */
1515: /*---------------------------------*/
1516: for (v=vStart; v<vEnd; v++) {
1517: /* Get row indices */
1518: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1519: DMNetworkGetNumVariables(dm,v,&nrows);
1520: if (!nrows) continue;
1522: DMNetworkIsGhostVertex(dm,v,&ghost);
1523: if (ghost) {
1524: PetscMalloc1(nrows,&rows_v);
1525: } else {
1526: rows_v = rows;
1527: }
1528: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
1530: /* Get supporting edges and connected vertices */
1531: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1533: for (e=0; e<nedges; e++) {
1534: /* Supporting edges */
1535: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1536: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1538: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1539: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1541: /* Connected vertices */
1542: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1543: vc = (v == cone[0]) ? cone[1]:cone[0];
1545: DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1546: DMNetworkGetNumVariables(dm,vc,&ncols);
1548: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1549: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1550: }
1552: /* Set matrix entries for vertex self */
1553: if (!ghost) {
1554: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1555: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1556: MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
1557: }
1558: if (ghost) {
1559: PetscFree(rows_v);
1560: }
1561: }
1562: PetscFree(rows);
1564: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1565: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1567: MatSetDM(*J,dm);
1568: return(0);
1569: }
1571: PetscErrorCode DMDestroy_Network(DM dm)
1572: {
1574: DM_Network *network = (DM_Network*) dm->data;
1577: if (--network->refct > 0) return(0);
1578: if (network->Je) {
1579: PetscFree(network->Je);
1580: }
1581: if (network->Jv) {
1582: PetscFree(network->Jvptr);
1583: PetscFree(network->Jv);
1584: }
1586: ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
1587: PetscSectionDestroy(&network->vertex.DofSection);
1588: PetscSectionDestroy(&network->vertex.GlobalDofSection);
1589: if (network->vertex.sf) {
1590: PetscSFDestroy(&network->vertex.sf);
1591: }
1592: /* edge */
1593: ISLocalToGlobalMappingDestroy(&network->edge.mapping);
1594: PetscSectionDestroy(&network->edge.DofSection);
1595: PetscSectionDestroy(&network->edge.GlobalDofSection);
1596: if (network->edge.sf) {
1597: PetscSFDestroy(&network->edge.sf);
1598: }
1599: DMDestroy(&network->plex);
1600: network->edges = NULL;
1601: PetscSectionDestroy(&network->DataSection);
1602: PetscSectionDestroy(&network->DofSection);
1604: PetscFree(network->componentdataarray);
1605: PetscFree(network->cvalue);
1606: PetscFree(network->header);
1607: PetscFree(network);
1608: return(0);
1609: }
1611: PetscErrorCode DMView_Network(DM dm, PetscViewer viewer)
1612: {
1614: DM_Network *network = (DM_Network*) dm->data;
1617: DMView(network->plex,viewer);
1618: return(0);
1619: }
1621: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
1622: {
1624: DM_Network *network = (DM_Network*) dm->data;
1627: DMGlobalToLocalBegin(network->plex,g,mode,l);
1628: return(0);
1629: }
1631: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
1632: {
1634: DM_Network *network = (DM_Network*) dm->data;
1637: DMGlobalToLocalEnd(network->plex,g,mode,l);
1638: return(0);
1639: }
1641: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
1642: {
1644: DM_Network *network = (DM_Network*) dm->data;
1647: DMLocalToGlobalBegin(network->plex,l,mode,g);
1648: return(0);
1649: }
1651: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
1652: {
1654: DM_Network *network = (DM_Network*) dm->data;
1657: DMLocalToGlobalEnd(network->plex,l,mode,g);
1658: return(0);
1659: }