Actual source code: network.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/
2: #include <petscdmplex.h>
3: #include <petscsf.h>
7: /*@
8: DMNetworkSetSizes - Sets the local and global vertices and edges.
10: Collective on DM
11:
12: Input Parameters:
13: + dm - the dm object
14: . nV - number of local vertices
15: . nE - number of local edges
16: . NV - number of global vertices (or PETSC_DETERMINE)
17: - NE - number of global edges (or PETSC_DETERMINE)
19: Notes
20: If one processor calls this with NV (NE) of PETSC_DECIDE then all processors must, otherwise the prgram will hang.
22: You cannot change the sizes once they have been set
24: Level: intermediate
26: .seealso: DMNetworkCreate
27: @*/
28: PetscErrorCode DMNetworkSetSizes(DM dm, PetscInt nV, PetscInt nE, PetscInt NV, PetscInt NE)
29: {
31: DM_Network *network = (DM_Network*) dm->data;
32: PetscInt a[2],b[2];
38: 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);
39: 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);
40: if ((network->nNodes >= 0 || network->NNodes >= 0) && (network->nNodes != nV || network->NNodes != 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->nNodes,network->NNodes);
41: 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);
42: if (NE < 0 || NV < 0) {
43: a[0] = nV; a[1] = nE;
44: MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
45: NV = b[0]; NE = b[1];
46: }
47: network->nNodes = nV;
48: network->NNodes = NV;
49: network->nEdges = nE;
50: network->NEdges = NE;
51: return(0);
52: }
56: /*@
57: DMNetworkSetEdgeList - Sets the list of local edges (vertex connectivity) for the network
59: Logically collective on DM
61: Input Parameters:
62: . edges - list of edges
64: Notes:
65: There is no copy involved in this operation, only the pointer is referenced. The edgelist should
66: not be destroyed before the call to DMNetworkLayoutSetUp
68: Level: intermediate
70: .seealso: DMNetworkCreate, DMNetworkSetSizes
71: @*/
72: PetscErrorCode DMNetworkSetEdgeList(DM dm, int edgelist[])
73: {
74: DM_Network *network = (DM_Network*) dm->data;
75:
77: network->edges = edgelist;
78: return(0);
79: }
83: /*@
84: DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network
86: Collective on DM
88: Input Parameters
89: . DM - the dmnetwork object
91: Notes:
92: This routine should be called after the network sizes and edgelists have been provided. It creates
93: the bare layout of the network and sets up the network to begin insertion of components.
95: All the components should be registered before calling this routine.
97: Level: intermediate
99: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
100: @*/
101: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
102: {
104: DM_Network *network = (DM_Network*) dm->data;
105: PetscInt dim = 1; /* One dimensional network */
106: PetscInt numCorners=2;
107: PetscInt spacedim=2;
108: double *vertexcoords=NULL;
109: PetscInt i;
110: PetscInt ndata;
113: if (network->nNodes) {
114: PetscCalloc1(numCorners*network->nNodes,&vertexcoords);
115: }
116: DMPlexCreateFromCellList(PetscObjectComm((PetscObject)dm),dim,network->nEdges,network->nNodes,numCorners,PETSC_FALSE,network->edges,spacedim,vertexcoords,&network->plex);
117: if (network->nNodes) {
118: PetscFree(vertexcoords);
119: }
120: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
121: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
122: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
124: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DataSection);
125: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DofSection);
126: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
127: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
129: network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
130: PetscCalloc1(network->pEnd-network->pStart,&network->header);
131: for (i = network->pStart; i < network->pEnd; i++) {
132: network->header[i].ndata = 0;
133: ndata = network->header[i].ndata;
134: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
135: network->header[i].offset[ndata] = 0;
136: }
137: PetscMalloc1(network->pEnd-network->pStart,&network->cvalue);
138: return(0);
139: }
143: /*@
144: DMNetworkRegisterComponent - Registers the network component
146: Logically collective on DM
148: Input Parameters
149: + dm - the network object
150: . name - the component name
151: - size - the storage size in bytes for this component data
153: Output Parameters
154: . key - an integer key that defines the component
156: Notes
157: This routine should be called by all processors before calling DMNetworkLayoutSetup().
159: Level: intermediate
161: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
162: @*/
163: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,PetscInt size,PetscInt *key)
164: {
165: PetscErrorCode ierr;
166: DM_Network *network = (DM_Network*) dm->data;
167: DMNetworkComponent *component=&network->component[network->ncomponent];
168: PetscBool flg=PETSC_FALSE;
169: PetscInt i;
173: for (i=0; i < network->ncomponent; i++) {
174: PetscStrcmp(component->name,name,&flg);
175: if (flg) {
176: *key = i;
177: return(0);
178: }
179: }
180:
181: PetscStrcpy(component->name,name);
182: component->size = size/sizeof(DMNetworkComponentGenericDataType);
183: *key = network->ncomponent;
184: network->ncomponent++;
185: return(0);
186: }
190: /*@
191: DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices.
193: Not Collective
195: Input Parameters:
196: + dm - The DMNetwork object
198: Output Paramters:
199: + vStart - The first vertex point
200: - vEnd - One beyond the last vertex point
202: Level: intermediate
204: .seealso: DMNetworkGetEdgeRange
205: @*/
206: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
207: {
208: DM_Network *network = (DM_Network*)dm->data;
211: if (vStart) *vStart = network->vStart;
212: if (vEnd) *vEnd = network->vEnd;
213: return(0);
214: }
218: /*@
219: DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges.
221: Not Collective
223: Input Parameters:
224: + dm - The DMNetwork object
226: Output Paramters:
227: + eStart - The first edge point
228: - eEnd - One beyond the last edge point
230: Level: intermediate
232: .seealso: DMNetworkGetVertexRange
233: @*/
234: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
235: {
236: DM_Network *network = (DM_Network*)dm->data;
239: if (eStart) *eStart = network->eStart;
240: if (eEnd) *eEnd = network->eEnd;
241: return(0);
242: }
246: /*@
247: DMNetworkAddComponent - Adds a network component at the given point (vertex/edge)
249: Not Collective
251: Input Parameters:
252: + dm - The DMNetwork object
253: . p - vertex/edge point
254: . componentkey - component key returned while registering the component
255: - compvalue - pointer to the data structure for the component
257: Level: intermediate
259: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
260: @*/
261: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
262: {
263: DM_Network *network = (DM_Network*)dm->data;
264: DMNetworkComponent *component = &network->component[componentkey];
265: DMNetworkComponentHeader header = &network->header[p];
266: DMNetworkComponentValue cvalue = &network->cvalue[p];
267: PetscErrorCode ierr;
268:
270: header->size[header->ndata] = component->size;
271: PetscSectionAddDof(network->DataSection,p,component->size);
272: header->key[header->ndata] = componentkey;
273: if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];
275: cvalue->data[header->ndata] = (void*)compvalue;
276: header->ndata++;
277: return(0);
278: }
282: /*@
283: DMNetworkGetNumComponents - Get the number of components at a vertex/edge
285: Not Collective
287: Input Parameters:
288: + dm - The DMNetwork object
289: . p - vertex/edge point
291: Output Parameters:
292: . numcomponents - Number of components at the vertex/edge
294: Level: intermediate
296: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
297: @*/
298: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
299: {
301: PetscInt offset;
302: DM_Network *network = (DM_Network*)dm->data;
305: PetscSectionGetOffset(network->DataSection,p,&offset);
306: *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
307: return(0);
308: }
312: /*@
313: DMNetworkGetComponentTypeOffset - Gets the type along with the offset for indexing the
314: component value from the component data array
316: Not Collective
318: Input Parameters:
319: + dm - The DMNetwork object
320: . p - vertex/edge point
321: - compnum - component number
322:
323: Output Parameters:
324: + compkey - the key obtained when registering the component
325: - offset - offset into the component data array associated with the vertex/edge point
327: Notes:
328: Typical usage:
330: DMNetworkGetComponentDataArray(dm, &arr);
331: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
332: Loop over vertices or edges
333: DMNetworkGetNumComponents(dm,v,&numcomps);
334: Loop over numcomps
335: DMNetworkGetComponentTypeOffset(dm,v,compnum,&key,&offset);
336: compdata = (UserCompDataType)(arr+offset);
337:
338: Level: intermediate
340: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
341: @*/
342: PetscErrorCode DMNetworkGetComponentTypeOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
343: {
344: PetscErrorCode ierr;
345: PetscInt offsetp;
346: DMNetworkComponentHeader header;
347: DM_Network *network = (DM_Network*)dm->data;
350: PetscSectionGetOffset(network->DataSection,p,&offsetp);
351: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
352: if (compkey) *compkey = header->key[compnum];
353: if (offset) *offset = offsetp+network->dataheadersize+header->offset[compnum];
354: return(0);
355: }
359: /*@
360: DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector.
362: Not Collective
364: Input Parameters:
365: + dm - The DMNetwork object
366: - p - the edge/vertex point
368: Output Parameters:
369: . offset - the offset
371: Level: intermediate
373: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
374: @*/
375: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
376: {
378: DM_Network *network = (DM_Network*)dm->data;
381: PetscSectionGetOffset(network->DofSection,p,offset);
382: return(0);
383: }
387: /*@
388: DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector.
390: Not Collective
392: Input Parameters:
393: + dm - The DMNetwork object
394: - p - the edge/vertex point
396: Output Parameters:
397: . offsetg - the offset
399: Level: intermediate
401: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
402: @*/
403: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
404: {
406: DM_Network *network = (DM_Network*)dm->data;
409: PetscSectionGetOffset(network->GlobalDofSection,p,offsetg);
410: if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost node */
411: return(0);
412: }
416: /*@
417: DMNetworkAddNumVariables - Add number of variables associated with a given point.
419: Not Collective
421: Input Parameters:
422: + dm - The DMNetworkObject
423: . p - the vertex/edge point
424: - nvar - number of additional variables
426: Level: intermediate
428: .seealso: DMNetworkSetNumVariables
429: @*/
430: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
431: {
433: DM_Network *network = (DM_Network*)dm->data;
436: PetscSectionAddDof(network->DofSection,p,nvar);
437: return(0);
438: }
442: /*@
443: DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point.
445: Not Collective
447: Input Parameters:
448: + dm - The DMNetworkObject
449: - p - the vertex/edge point
451: Output Parameters:
452: . nvar - number of variables
454: Level: intermediate
456: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
457: @*/
458: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
459: {
461: DM_Network *network = (DM_Network*)dm->data;
464: PetscSectionGetDof(network->DofSection,p,nvar);
465: return(0);
466: }
470: /*@
471: DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point.
473: Not Collective
475: Input Parameters:
476: + dm - The DMNetworkObject
477: . p - the vertex/edge point
478: - nvar - number of variables
480: Level: intermediate
482: .seealso: DMNetworkAddNumVariables
483: @*/
484: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
485: {
487: DM_Network *network = (DM_Network*)dm->data;
490: PetscSectionSetDof(network->DofSection,p,nvar);
491: return(0);
492: }
494: /* Sets up the array that holds the data for all components and its associated section. This
495: function is called during DMSetUp() */
498: PetscErrorCode DMNetworkComponentSetUp(DM dm)
499: {
500: PetscErrorCode ierr;
501: DM_Network *network = (DM_Network*)dm->data;
502: PetscInt arr_size;
503: PetscInt p,offset,offsetp;
504: DMNetworkComponentHeader header;
505: DMNetworkComponentValue cvalue;
506: DMNetworkComponentGenericDataType *componentdataarray;
507: PetscInt ncomp, i;
510: PetscSectionSetUp(network->DataSection);
511: PetscSectionGetStorageSize(network->DataSection,&arr_size);
512: PetscMalloc1(arr_size,&network->componentdataarray);
513: componentdataarray = network->componentdataarray;
514: for (p = network->pStart; p < network->pEnd; p++) {
515: PetscSectionGetOffset(network->DataSection,p,&offsetp);
516: /* Copy header */
517: header = &network->header[p];
518: PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
519: /* Copy data */
520: cvalue = &network->cvalue[p];
521: ncomp = header->ndata;
522: for (i = 0; i < ncomp; i++) {
523: offset = offsetp + network->dataheadersize + header->offset[i];
524: PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
525: }
526: }
527: return(0);
528: }
530: /* Sets up the section for dofs. This routine is called during DMSetUp() */
533: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
534: {
536: DM_Network *network = (DM_Network*)dm->data;
539: PetscSectionSetUp(network->DofSection);
540: return(0);
541: }
545: /*@C
546: DMNetworkGetComponentDataArray - Returns the component data array
548: Not Collective
550: Input Parameters:
551: . dm - The DMNetwork Object
553: Output Parameters:
554: . componentdataarray - array that holds data for all components
556: Level: intermediate
558: .seealso: DMNetworkGetComponentTypeOffset, DMNetworkGetNumComponents
559: @*/
560: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
561: {
562: DM_Network *network = (DM_Network*)dm->data;
565: *componentdataarray = network->componentdataarray;
566: return(0);
567: }
571: /*@
572: DMNetworkDistribute - Distributes the network and moves associated component data.
574: Collective
576: Input Parameter:
577: + oldDM - the original DMNetwork object
578: - overlap - The overlap of partitions, 0 is the default
580: Output Parameter:
581: . distDM - the distributed DMNetwork object
583: Notes:
584: This routine should be called only when using multiple processors.
586: Distributes the network with <overlap>-overlapping partitioning of the edges.
588: Level: intermediate
590: .seealso: DMNetworkCreate
591: @*/
592: PetscErrorCode DMNetworkDistribute(DM oldDM, PetscInt overlap,DM *distDM)
593: {
595: DM_Network *oldDMnetwork = (DM_Network*)oldDM->data;
596: PetscSF pointsf;
597: DM newDM;
598: DM_Network *newDMnetwork;
601: DMNetworkCreate(PetscObjectComm((PetscObject)oldDM),&newDM);
602: newDMnetwork = (DM_Network*)newDM->data;
603: newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
604: /* Distribute plex dm and dof section */
605: DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
606: /* Distribute dof section */
607: PetscSectionCreate(PetscObjectComm((PetscObject)oldDM),&newDMnetwork->DofSection);
608: PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
609: PetscSectionCreate(PetscObjectComm((PetscObject)oldDM),&newDMnetwork->DataSection);
610: /* Distribute data and associated section */
611: DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
612: /* Destroy point SF */
613: PetscSFDestroy(&pointsf);
614:
615: PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
616: DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
617: DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
618: newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
619: newDMnetwork->nNodes = newDMnetwork->vEnd - newDMnetwork->vStart;
620: newDMnetwork->NNodes = oldDMnetwork->NNodes;
621: newDMnetwork->NEdges = oldDMnetwork->NEdges;
622: /* Set Dof section as the default section for dm */
623: DMSetDefaultSection(newDMnetwork->plex,newDMnetwork->DofSection);
624: DMGetDefaultGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);
626: *distDM = newDM;
627: return(0);
628: }
632: /*@C
633: DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point
635: Not Collective
637: Input Parameters:
638: + dm - The DMNetwork object
639: - p - the vertex point
641: Output Paramters:
642: + nedges - number of edges connected to this vertex point
643: - edges - List of edge points
645: Level: intermediate
647: Fortran Notes:
648: Since it returns an array, this routine is only available in Fortran 90, and you must
649: include petsc.h90 in your code.
651: .seealso: DMNetworkCreate, DMNetworkGetConnectedNodes
652: @*/
653: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
654: {
656: DM_Network *network = (DM_Network*)dm->data;
659: DMPlexGetSupportSize(network->plex,vertex,nedges);
660: DMPlexGetSupport(network->plex,vertex,edges);
661: return(0);
662: }
666: /*@C
667: DMNetworkGetConnectedNodes - Return the connected vertices for this edge point
669: Not Collective
671: Input Parameters:
672: + dm - The DMNetwork object
673: - p - the edge point
675: Output Paramters:
676: . vertices - vertices connected to this edge
678: Level: intermediate
680: Fortran Notes:
681: Since it returns an array, this routine is only available in Fortran 90, and you must
682: include petsc.h90 in your code.
684: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
685: @*/
686: PetscErrorCode DMNetworkGetConnectedNodes(DM dm,PetscInt edge,const PetscInt *vertices[])
687: {
689: DM_Network *network = (DM_Network*)dm->data;
692: DMPlexGetCone(network->plex,edge,vertices);
693: return(0);
694: }
698: /*@
699: DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex
701: Not Collective
703: Input Parameters:
704: + dm - The DMNetwork object
705: . p - the vertex point
707: Output Parameter:
708: . isghost - TRUE if the vertex is a ghost point
710: Level: intermediate
712: .seealso: DMNetworkCreate, DMNetworkGetConnectedNodes, DMNetworkGetVertexRange
713: @*/
714: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
715: {
717: DM_Network *network = (DM_Network*)dm->data;
718: PetscInt offsetg;
719: PetscSection sectiong;
722: *isghost = PETSC_FALSE;
723: DMGetDefaultGlobalSection(network->plex,§iong);
724: PetscSectionGetOffset(sectiong,p,&offsetg);
725: if (offsetg < 0) *isghost = PETSC_TRUE;
726: return(0);
727: }
731: PetscErrorCode DMSetUp_Network(DM dm)
732: {
734: DM_Network *network=(DM_Network*)dm->data;
737: DMNetworkComponentSetUp(dm);
738: DMNetworkVariablesSetUp(dm);
740: DMSetDefaultSection(network->plex,network->DofSection);
741: DMGetDefaultGlobalSection(network->plex,&network->GlobalDofSection);
742: return(0);
743: }
747: /*@
748: DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
749: -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?
751: Collective
753: Input Parameters:
754: + dm - The DMNetwork object
755: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
756: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices
758: Level: intermediate
760: @*/
761: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
762: {
763: DM_Network *network=(DM_Network*)dm->data;
766: network->userEdgeJacobian = eflg;
767: network->userVertexJacobian = vflg;
768: return(0);
769: }
773: /*@
774: DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network
776: Not Collective
778: Input Parameters:
779: + dm - The DMNetwork object
780: . p - the edge point
781: - J - array (size = 3) of Jacobian submatrices for this edge point:
782: J[0]: this edge
783: J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedNodes()
785: Level: intermediate
787: .seealso: DMNetworkVertexSetMatrix
788: @*/
789: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
790: {
792: DM_Network *network=(DM_Network*)dm->data;
795: if (!network->userEdgeJacobian) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
796: if (!network->Je) {
797: PetscCalloc1(3*network->nEdges,&network->Je);
798: }
799: network->Je[3*p] = J[0];
800: network->Je[3*p+1] = J[1];
801: network->Je[3*p+2] = J[2];
802: return(0);
803: }
807: /*@
808: DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network
810: Not Collective
812: Input Parameters:
813: + dm - The DMNetwork object
814: . p - the vertex point
815: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
816: J[0]: this vertex
817: J[1+2*i]: i-th supporting edge
818: J[1+2*i+1]: i-th connected vertex
820: Level: intermediate
822: .seealso: DMNetworkEdgeSetMatrix
823: @*/
824: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
825: {
827: DM_Network *network=(DM_Network*)dm->data;
828: PetscInt i,*vptr,nedges,vStart,vEnd;
829: const PetscInt *edges;
832: if (!network->userVertexJacobian) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");
834: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
836: if (!network->Jv) {
837: PetscInt nNodes = network->nNodes,nedges_total;
838: /* count nvertex_total */
839: nedges_total = 0;
840: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
841: PetscMalloc1(nNodes+1,&vptr);
842:
843: vptr[0] = 0;
844: for (i=0; i<nNodes; i++) {
845: DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
846: nedges_total += nedges;
847: vptr[i+1] = vptr[i] + 2*nedges + 1;
848: }
850: PetscCalloc1(2*nedges_total+nNodes,&network->Jv);
851: network->Jvptr = vptr;
852: }
854: vptr = network->Jvptr;
855: network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */
857: /* Set Jacobian for each supporting edge and connected vertex */
858: DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
859: for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
860: return(0);
861: }
865: PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
866: {
868: PetscInt j,*cols;
869: PetscScalar *zeros;
872: PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
873: for (j=0; j<ncols; j++) cols[j] = j+ cstart;
874: MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
875: PetscFree2(cols,zeros);
876: return(0);
877: }
881: PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
882: {
884: PetscInt j,M,N,row,col,ncols_u;
885: const PetscInt *cols;
886: PetscScalar zero=0.0;
889: MatGetSize(Ju,&M,&N);
890: if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);
892: for (row=0; row<nrows; row++) {
893: MatGetRow(Ju,row,&ncols_u,&cols,NULL);
894: for (j=0; j<ncols_u; j++) {
895: col = cols[j] + cstart;
896: MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
897: }
898: MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
899: }
900: return(0);
901: }
905: PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
906: {
910: if (Ju) {
911: MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
912: } else {
913: MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
914: }
915: return(0);
916: }
920: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
921: {
923: DM_Network *network = (DM_Network*) dm->data;
924: PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
925: PetscInt cstart,ncols,j,e,v,*dnz,*onz,*dnzu,*onzu;
926: PetscBool ghost;
927: Mat Juser;
928: PetscSection sectionGlobal;
929: PetscInt nedges,*vptr=NULL,vc; /* suppress maybe-uninitialized warning */
930: const PetscInt *edges,*cone;
933: if (!network->userEdgeJacobian && !network->userVertexJacobian) {
934: /* user does not provide Jacobian blocks */
935: DMCreateMatrix(network->plex,J);
936: MatSetDM(*J,dm);
937: return(0);
938: }
939:
940: MatCreate(PetscObjectComm((PetscObject)dm),J);
941: DMGetDefaultGlobalSection(network->plex,§ionGlobal);
942: PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
943: MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
945: MatSetType(*J,MATAIJ);
946: MatSetFromOptions(*J);
948: /* Preallocation - submatrix for an element (edge/vertex) is allocated as a dense block,
949: see DMCreateMatrix_Plex() */
950: PetscCalloc4(localSize,&dnz,localSize,&onz,localSize,&dnzu,localSize,&onzu);
951: DMPlexPreallocateOperator(network->plex,1,dnz,onz,dnzu,onzu,*J,PETSC_FALSE);
952: PetscFree4(dnz,onz,dnzu,onzu);
954: /* Set matrix entries for edges */
955: /*------------------------------*/
956: DMNetworkGetEdgeRange(dm,&eStart,&eEnd);
958: for (e=eStart; e<eEnd; e++) {
959: /* Get row indices */
960: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
961: DMNetworkGetNumVariables(dm,e,&nrows);
962: if (nrows) {
963: if (!network->Je) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Je");
964: PetscMalloc1(nrows,&rows);
965: for (j=0; j<nrows; j++) rows[j] = j + rstart;
967: /* Set matrix entries for conntected vertices */
968: DMNetworkGetConnectedNodes(dm,e,&cone);
969: for (v=0; v<2; v++) {
970: DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
971: DMNetworkGetNumVariables(dm,cone[v],&ncols);
973: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
974: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
975: }
977: /* Set matrix entries for edge self */
978: cstart = rstart;
979: Juser = network->Je[3*e]; /* Jacobian(e,e) */
980: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
982: PetscFree(rows);
983: }
984: }
986: /* Set matrix entries for vertices */
987: /*---------------------------------*/
988: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
989: if (vEnd - vStart) {
990: if (!network->Jv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Jv");
991: vptr = network->Jvptr;
992: if (!vptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide vptr");
993: }
994:
995: for (v=vStart; v<vEnd; v++) {
996: /* Get row indices */
997: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
998: DMNetworkGetNumVariables(dm,v,&nrows);
999: if (!nrows) continue;
1001: PetscMalloc1(nrows,&rows);
1002: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1004: /* Get supporting edges and connected vertices */
1005: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1006:
1007: for (e=0; e<nedges; e++) {
1008: /* Supporting edges */
1009: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1010: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1012: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1013: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
1015: /* Connected vertices */
1016: DMNetworkGetConnectedNodes(dm,edges[e],&cone);
1017: vc = (v == cone[0]) ? cone[1]:cone[0];
1018:
1019: DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1020: DMNetworkGetNumVariables(dm,vc,&ncols);
1022: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1023: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1024: }
1026: /* Set matrix entries for vertex self */
1027: DMNetworkIsGhostVertex(dm,v,&ghost);
1028: if (!ghost) {
1029: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1030: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1031: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1032: }
1033: PetscFree(rows);
1034: }
1035: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1036: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1038: MatSetDM(*J,dm);
1039: return(0);
1040: }
1044: PetscErrorCode DMDestroy_Network(DM dm)
1045: {
1047: DM_Network *network = (DM_Network*) dm->data;
1050: if (--network->refct > 0) return(0);
1051: if (network->Je) {
1052: PetscFree(network->Je);
1053: }
1054: if (network->Jv) {
1055: PetscFree(network->Jvptr);
1056: PetscFree(network->Jv);
1057: }
1058: DMDestroy(&network->plex);
1059: network->edges = NULL;
1060: PetscSectionDestroy(&network->DataSection);
1061: PetscSectionDestroy(&network->DofSection);
1062:
1063: PetscFree(network->componentdataarray);
1064: PetscFree(network->cvalue);
1065: PetscFree(network->header);
1066: PetscFree(network);
1067: return(0);
1068: }
1072: PetscErrorCode DMView_Network(DM dm, PetscViewer viewer)
1073: {
1075: DM_Network *network = (DM_Network*) dm->data;
1078: DMView(network->plex,viewer);
1079: return(0);
1080: }
1084: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
1085: {
1087: DM_Network *network = (DM_Network*) dm->data;
1090: DMGlobalToLocalBegin(network->plex,g,mode,l);
1091: return(0);
1092: }
1096: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
1097: {
1099: DM_Network *network = (DM_Network*) dm->data;
1102: DMGlobalToLocalEnd(network->plex,g,mode,l);
1103: return(0);
1104: }
1108: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
1109: {
1111: DM_Network *network = (DM_Network*) dm->data;
1114: DMLocalToGlobalBegin(network->plex,l,mode,g);
1115: return(0);
1116: }
1120: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
1121: {
1123: DM_Network *network = (DM_Network*) dm->data;
1126: DMLocalToGlobalEnd(network->plex,l,mode,g);
1127: return(0);
1128: }