Actual source code: dagetelem.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
6: static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
7: {
9: DM_DA *da = (DM_DA*)dm->data;
10: PetscInt i,xs,xe,Xs,Xe;
11: PetscInt cnt=0;
14: if (!da->e) {
15: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
16: DMDAGetCorners(dm,&xs,0,0,&xe,0,0);
17: DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);
18: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
19: da->ne = 1*(xe - xs - 1);
20: PetscMalloc1(1 + 2*da->ne,&da->e);
21: for (i=xs; i<xe-1; i++) {
22: da->e[cnt++] = (i-Xs);
23: da->e[cnt++] = (i-Xs+1);
24: }
25: }
26: *nel = da->ne;
27: *nen = 2;
28: *e = da->e;
29: return(0);
30: }
34: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
35: {
37: DM_DA *da = (DM_DA*)dm->data;
38: PetscInt i,xs,xe,Xs,Xe;
39: PetscInt j,ys,ye,Ys,Ye;
40: PetscInt cnt=0, cell[4], ns=2, nn=3;
41: PetscInt c, split[] = {0,1,3,
42: 2,3,1};
45: if (!da->e) {
46: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
47: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;}
48: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;}
49: DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
50: DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
51: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
52: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
53: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
54: PetscMalloc1(1 + nn*da->ne,&da->e);
55: for (j=ys; j<ye-1; j++) {
56: for (i=xs; i<xe-1; i++) {
57: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs);
58: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs);
59: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
60: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs);
61: if (da->elementtype == DMDA_ELEMENT_P1) {
62: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
63: }
64: if (da->elementtype == DMDA_ELEMENT_Q1) {
65: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
66: }
67: }
68: }
69: }
70: *nel = da->ne;
71: *nen = nn;
72: *e = da->e;
73: return(0);
74: }
78: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
79: {
81: DM_DA *da = (DM_DA*)dm->data;
82: PetscInt i,xs,xe,Xs,Xe;
83: PetscInt j,ys,ye,Ys,Ye;
84: PetscInt k,zs,ze,Zs,Ze;
85: PetscInt cnt=0, cell[8], ns=6, nn=4;
86: PetscInt c, split[] = {0,1,3,7,
87: 0,1,7,4,
88: 1,2,3,7,
89: 1,2,7,6,
90: 1,4,5,7,
91: 1,5,6,7};
94: if (!da->e) {
95: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
96: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;}
97: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;}
98: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
99: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
100: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
101: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
102: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
103: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
104: PetscMalloc1(1 + nn*da->ne,&da->e);
105: for (k=zs; k<ze-1; k++) {
106: for (j=ys; j<ye-1; j++) {
107: for (i=xs; i<xe-1; i++) {
108: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
109: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
110: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
111: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
112: cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
113: cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
114: cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
115: cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
116: if (da->elementtype == DMDA_ELEMENT_P1) {
117: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
118: }
119: if (da->elementtype == DMDA_ELEMENT_Q1) {
120: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
121: }
122: }
123: }
124: }
125: }
126: *nel = da->ne;
127: *nen = nn;
128: *e = da->e;
129: return(0);
130: }
132: /*@C
133: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
135: Not Collective
137: Input Parameter:
138: . da - the DMDA object
140: Output Parameters:
141: . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1
143: Level: intermediate
145: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
146: @*/
149: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
150: {
151: DM_DA *dd = (DM_DA*)da->data;
157: if (dd->elementtype != etype) {
158: PetscFree(dd->e);
160: dd->elementtype = etype;
161: dd->ne = 0;
162: dd->e = NULL;
163: }
164: return(0);
165: }
167: /*@C
168: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
170: Not Collective
172: Input Parameter:
173: . da - the DMDA object
175: Output Parameters:
176: . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1
178: Level: intermediate
180: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
181: @*/
184: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
185: {
186: DM_DA *dd = (DM_DA*)da->data;
191: *etype = dd->elementtype;
192: return(0);
193: }
195: /*@C
196: DMDAGetElements - Gets an array containing the indices (in local coordinates)
197: of all the local elements
199: Not Collective
201: Input Parameter:
202: . dm - the DM object
204: Output Parameters:
205: + nel - number of local elements
206: . nen - number of element nodes
207: - e - the local indices of the elements' vertices
209: Level: intermediate
211: Notes:
212: Call DMDARestoreElements() once you have finished accessing the elements.
214: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
216: If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.
218: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
219: @*/
222: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
223: {
224: PetscInt dim;
226: DM_DA *dd = (DM_DA*)dm->data;
229: if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX");
230: DMGetDimension(dm, &dim);
231: if (dim==-1) {
232: *nel = 0; *nen = 0; *e = NULL;
233: } else if (dim==1) {
234: DMDAGetElements_1D(dm,nel,nen,e);
235: } else if (dim==2) {
236: DMDAGetElements_2D(dm,nel,nen,e);
237: } else if (dim==3) {
238: DMDAGetElements_3D(dm,nel,nen,e);
239: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
240: return(0);
241: }
245: /*@C
246: DMDARestoreElements - Restores the array obtained with DMDAGetElements()
248: Not Collective
250: Input Parameter:
251: + dm - the DM object
252: . nel - number of local elements
253: . nen - number of element nodes
254: - e - the local indices of the elements' vertices
256: Level: intermediate
258: Note: You should not access these values after you have called this routine.
260: This restore signals the DMDA object that you no longer need access to the array information.
262: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
263: @*/
264: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
265: {
271: *nel = 0;
272: *nen = -1;
273: *e = NULL;
274: return(0);
275: }