Actual source code: dagetelem.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/daimpl.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;
13: if (!da->e) {
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: PetscMalloc((1 + 2*da->ne)*sizeof(PetscInt),&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: }
32: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
33: {
35: DM_DA *da = (DM_DA*)dm->data;
36: PetscInt i,xs,xe,Xs,Xe;
37: PetscInt j,ys,ye,Ys,Ye;
38: PetscInt cnt=0, cell[4], ns=2, nn=3;
39: PetscInt c, split[] = {0,1,3,
40: 2,3,1};
42: if (!da->e) {
43: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;}
44: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;}
45: DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
46: DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
47: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
48: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
49: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
50: PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);
51: for (j=ys; j<ye-1; j++) {
52: for (i=xs; i<xe-1; i++) {
53: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs);
54: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs);
55: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
56: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs);
57: if (da->elementtype == DMDA_ELEMENT_P1) {
58: for (c=0; c<ns*nn; c++)
59: da->e[cnt++] = cell[split[c]];
60: }
61: if (da->elementtype == DMDA_ELEMENT_Q1) {
62: for (c=0; c<ns*nn; c++)
63: 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};
91: if (!da->e) {
92: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;}
93: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;}
94: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
95: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
96: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
97: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
98: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
99: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
100: PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);
101: for (k=zs; k<ze-1; k++) {
102: for (j=ys; j<ye-1; j++) {
103: for (i=xs; i<xe-1; i++) {
104: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
105: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
106: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
107: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
108: cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
109: cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
110: cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
111: cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
112: if (da->elementtype == DMDA_ELEMENT_P1) {
113: for (c=0; c<ns*nn; c++)
114: da->e[cnt++] = cell[split[c]];
115: }
116: if (da->elementtype == DMDA_ELEMENT_Q1) {
117: for (c=0; c<ns*nn; c++)
118: 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 ELEMENT_Q1
141: Level: intermediate
143: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
144: @*/
147: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
148: {
149: DM_DA *dd = (DM_DA*)da->data;
155: if (dd->elementtype != etype) {
156: PetscFree(dd->e);
157: dd->elementtype = etype;
158: dd->ne = 0;
159: dd->e = PETSC_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;
187: *etype = dd->elementtype;
188: return(0);
189: }
191: /*@C
192: DMDAGetElements - Gets an array containing the indices (in local coordinates)
193: of all the local elements
195: Not Collective
197: Input Parameter:
198: . dm - the DM object
200: Output Parameters:
201: + nel - number of local elements
202: . nen - number of element nodes
203: - e - the local indices of the elements' vertices
205: Level: intermediate
207: .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
208: @*/
211: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
212: {
213: DM_DA *da = (DM_DA*)dm->data;
216: if (da->dim==-1) {
217: *nel = 0; *nen = 0; *e = PETSC_NULL;
218: } else if (da->dim==1) {
219: DMDAGetElements_1D(dm,nel,nen,e);
220: } else if (da->dim==2) {
221: DMDAGetElements_2D(dm,nel,nen,e);
222: } else if (da->dim==3) {
223: DMDAGetElements_3D(dm,nel,nen,e);
224: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim);
226: return(0);
227: }
231: /*@C
232: DMDARestoreElements - Returns an array containing the indices (in local coordinates)
233: of all the local elements obtained with DMDAGetElements()
235: Not Collective
237: Input Parameter:
238: + dm - the DM object
239: . nel - number of local elements
240: . nen - number of element nodes
241: - e - the local indices of the elements' vertices
243: Level: intermediate
245: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
246: @*/
247: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
248: {
254: /* XXX */
255: return(0);
256: }