Actual source code: dagetelem.c
petsc-3.8.4 2018-03-24
2: #include <petsc/private/dmdaimpl.h>
4: static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
5: {
7: DM_DA *da = (DM_DA*)dm->data;
8: PetscInt i,xs,xe,Xs,Xe;
9: PetscInt cnt=0;
12: if (!da->e) {
13: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
14: DMDAGetCorners(dm,&xs,0,0,&xe,0,0);
15: DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);
16: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
17: da->ne = 1*(xe - xs - 1);
18: PetscMalloc1(1 + 2*da->ne,&da->e);
19: for (i=xs; i<xe-1; i++) {
20: da->e[cnt++] = (i-Xs);
21: da->e[cnt++] = (i-Xs+1);
22: }
23: }
24: *nel = da->ne;
25: *nen = 2;
26: *e = da->e;
27: return(0);
28: }
30: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
31: {
33: DM_DA *da = (DM_DA*)dm->data;
34: PetscInt i,xs,xe,Xs,Xe;
35: PetscInt j,ys,ye,Ys,Ye;
36: PetscInt cnt=0, cell[4], ns=2, nn=3;
37: PetscInt c, split[] = {0,1,3,
38: 2,3,1};
41: if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;}
42: if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;}
43: if (!da->e) {
44: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
45: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
46: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
47: DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
48: DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
49: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
50: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
51: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
52: PetscMalloc1(1 + nn*da->ne,&da->e);
53: for (j=ys; j<ye-1; j++) {
54: for (i=xs; i<xe-1; i++) {
55: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs);
56: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs);
57: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
58: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs);
59: if (da->elementtype == DMDA_ELEMENT_P1) {
60: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
61: }
62: if (da->elementtype == DMDA_ELEMENT_Q1) {
63: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
64: }
65: }
66: }
67: }
68: *nel = da->ne;
69: *nen = nn;
70: *e = da->e;
71: return(0);
72: }
74: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
75: {
77: DM_DA *da = (DM_DA*)dm->data;
78: PetscInt i,xs,xe,Xs,Xe;
79: PetscInt j,ys,ye,Ys,Ye;
80: PetscInt k,zs,ze,Zs,Ze;
81: PetscInt cnt=0, cell[8], ns=6, nn=4;
82: PetscInt c, split[] = {0,1,3,7,
83: 0,1,7,4,
84: 1,2,3,7,
85: 1,2,7,6,
86: 1,4,5,7,
87: 1,5,6,7};
90: if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;}
91: if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;}
92: if (!da->e) {
93: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
94: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
95: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
96: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
97: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
98: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
99: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
100: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
101: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
102: PetscMalloc1(1 + nn*da->ne,&da->e);
103: for (k=zs; k<ze-1; k++) {
104: for (j=ys; j<ye-1; j++) {
105: for (i=xs; i<xe-1; i++) {
106: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
107: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
108: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
109: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
110: cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
111: cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
112: cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
113: cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
114: if (da->elementtype == DMDA_ELEMENT_P1) {
115: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
116: }
117: if (da->elementtype == DMDA_ELEMENT_Q1) {
118: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
119: }
120: }
121: }
122: }
123: }
124: *nel = da->ne;
125: *nen = nn;
126: *e = da->e;
127: return(0);
128: }
130: /*@C
131: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
133: Not Collective
135: Input Parameter:
136: . da - the DMDA object
138: Output Parameters:
139: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
141: Level: intermediate
143: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
144: @*/
145: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
146: {
147: DM_DA *dd = (DM_DA*)da->data;
153: if (dd->elementtype != etype) {
154: PetscFree(dd->e);
156: dd->elementtype = etype;
157: dd->ne = 0;
158: dd->e = NULL;
159: }
160: return(0);
161: }
163: /*@C
164: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
166: Not Collective
168: Input Parameter:
169: . da - the DMDA object
171: Output Parameters:
172: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
174: Level: intermediate
176: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
177: @*/
178: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
179: {
180: DM_DA *dd = (DM_DA*)da->data;
185: *etype = dd->elementtype;
186: return(0);
187: }
189: /*@C
190: DMDAGetElements - Gets an array containing the indices (in local coordinates)
191: of all the local elements
193: Not Collective
195: Input Parameter:
196: . dm - the DM object
198: Output Parameters:
199: + nel - number of local elements
200: . nen - number of element nodes
201: - e - the local indices of the elements' vertices
203: Level: intermediate
205: Notes:
206: Call DMDARestoreElements() once you have finished accessing the elements.
208: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
210: If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.
212: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
213: @*/
214: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
215: {
216: PetscInt dim;
218: DM_DA *dd = (DM_DA*)dm->data;
221: 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");
222: DMGetDimension(dm, &dim);
223: if (dim==-1) {
224: *nel = 0; *nen = 0; *e = NULL;
225: } else if (dim==1) {
226: DMDAGetElements_1D(dm,nel,nen,e);
227: } else if (dim==2) {
228: DMDAGetElements_2D(dm,nel,nen,e);
229: } else if (dim==3) {
230: DMDAGetElements_3D(dm,nel,nen,e);
231: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
232: return(0);
233: }
235: /*@C
236: DMDARestoreElements - Restores the array obtained with DMDAGetElements()
238: Not Collective
240: Input Parameter:
241: + dm - the DM object
242: . nel - number of local elements
243: . nen - number of element nodes
244: - e - the local indices of the elements' vertices
246: Level: intermediate
248: Note: You should not access these values after you have called this routine.
250: This restore signals the DMDA object that you no longer need access to the array information.
252: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
253: @*/
254: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
255: {
261: *nel = 0;
262: *nen = -1;
263: *e = NULL;
264: return(0);
265: }