Actual source code: dagetelem.c
petsc-3.5.4 2015-05-23
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: DMDAGetCorners(dm,&xs,0,0,&xe,0,0);
16: DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);
17: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
18: da->ne = 1*(xe - xs - 1);
19: PetscMalloc1((1 + 2*da->ne),&da->e);
20: for (i=xs; i<xe-1; i++) {
21: da->e[cnt++] = (i-Xs);
22: da->e[cnt++] = (i-Xs+1);
23: }
24: }
25: *nel = da->ne;
26: *nen = 2;
27: *e = da->e;
28: return(0);
29: }
33: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
34: {
36: DM_DA *da = (DM_DA*)dm->data;
37: PetscInt i,xs,xe,Xs,Xe;
38: PetscInt j,ys,ye,Ys,Ye;
39: PetscInt cnt=0, cell[4], ns=2, nn=3;
40: PetscInt c, split[] = {0,1,3,
41: 2,3,1};
44: if (!da->e) {
45: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;}
46: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;}
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: }
76: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
77: {
79: DM_DA *da = (DM_DA*)dm->data;
80: PetscInt i,xs,xe,Xs,Xe;
81: PetscInt j,ys,ye,Ys,Ye;
82: PetscInt k,zs,ze,Zs,Ze;
83: PetscInt cnt=0, cell[8], ns=6, nn=4;
84: PetscInt c, split[] = {0,1,3,7,
85: 0,1,7,4,
86: 1,2,3,7,
87: 1,2,7,6,
88: 1,4,5,7,
89: 1,5,6,7};
92: if (!da->e) {
93: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;}
94: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;}
95: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
96: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
97: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
98: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
99: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
100: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
101: PetscMalloc1((1 + nn*da->ne),&da->e);
102: for (k=zs; k<ze-1; k++) {
103: for (j=ys; j<ye-1; j++) {
104: for (i=xs; i<xe-1; i++) {
105: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
106: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
107: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
108: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
109: cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
110: cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
111: cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
112: cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
113: if (da->elementtype == DMDA_ELEMENT_P1) {
114: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
115: }
116: if (da->elementtype == DMDA_ELEMENT_Q1) {
117: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
118: }
119: }
120: }
121: }
122: }
123: *nel = da->ne;
124: *nen = nn;
125: *e = da->e;
126: return(0);
127: }
129: /*@C
130: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
132: Not Collective
134: Input Parameter:
135: . da - the DMDA object
137: Output Parameters:
138: . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1
140: Level: intermediate
142: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
143: @*/
146: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
147: {
148: DM_DA *dd = (DM_DA*)da->data;
154: if (dd->elementtype != etype) {
155: PetscFree(dd->e);
157: dd->elementtype = etype;
158: dd->ne = 0;
159: dd->e = NULL;
160: }
161: return(0);
162: }
164: /*@C
165: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
167: Not Collective
169: Input Parameter:
170: . da - the DMDA object
172: Output Parameters:
173: . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1
175: Level: intermediate
177: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
178: @*/
181: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
182: {
183: DM_DA *dd = (DM_DA*)da->data;
188: *etype = dd->elementtype;
189: return(0);
190: }
192: /*@C
193: DMDAGetElements - Gets an array containing the indices (in local coordinates)
194: of all the local elements
196: Not Collective
198: Input Parameter:
199: . dm - the DM object
201: Output Parameters:
202: + nel - number of local elements
203: . nen - number of element nodes
204: - e - the local indices of the elements' vertices
206: Level: intermediate
208: .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
209: @*/
212: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
213: {
214: DM_DA *da = (DM_DA*)dm->data;
218: if (da->dim==-1) {
219: *nel = 0; *nen = 0; *e = NULL;
220: } else if (da->dim==1) {
221: DMDAGetElements_1D(dm,nel,nen,e);
222: } else if (da->dim==2) {
223: DMDAGetElements_2D(dm,nel,nen,e);
224: } else if (da->dim==3) {
225: DMDAGetElements_3D(dm,nel,nen,e);
226: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim);
227: return(0);
228: }
232: /*@C
233: DMDARestoreElements - Returns an array containing the indices (in local coordinates)
234: of all the local elements obtained with DMDAGetElements()
236: Not Collective
238: Input Parameter:
239: + dm - the DM object
240: . nel - number of local elements
241: . nen - number of element nodes
242: - e - the local indices of the elements' vertices
244: Level: intermediate
246: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
247: @*/
248: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
249: {
255: /* XXX */
256: return(0);
257: }