1 2 #include <private/daimpl.h> /*I "petscdmda.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMDAGetElements_1D" 6 static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 7 { 8 PetscErrorCode ierr; 9 DM_DA *da = (DM_DA*)dm->data; 10 PetscInt i,xs,xe,Xs,Xe; 11 PetscInt cnt=0; 12 PetscFunctionBegin; 13 if (!da->e) { 14 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 15 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 16 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 17 da->ne = 1*(xe - xs - 1); 18 ierr = PetscMalloc((1 + 2*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "DMDAGetElements_2D" 32 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 33 { 34 PetscErrorCode ierr; 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}; 41 PetscFunctionBegin; 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 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 46 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 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 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "DMDAGetElements_3D" 76 static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 77 { 78 PetscErrorCode ierr; 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}; 90 PetscFunctionBegin; 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 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 95 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 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 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 128 } 129 130 /*@C 131 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 132 133 Not Collective 134 135 Input Parameter: 136 . da - the DMDA object 137 138 Output Parameters: 139 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 140 141 Level: intermediate 142 143 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 144 @*/ 145 #undef __FUNCT__ 146 #define __FUNCT__ "DMDASetElementType" 147 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 148 { 149 DM_DA *dd = (DM_DA*)da->data; 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 PetscValidHeaderSpecific(da,DM_CLASSID,1); 154 PetscValidLogicalCollectiveEnum(da,etype,2); 155 if (dd->elementtype != etype) { 156 ierr = PetscFree(dd->e);CHKERRQ(ierr); 157 dd->elementtype = etype; 158 dd->ne = 0; 159 dd->e = PETSC_NULL; 160 } 161 PetscFunctionReturn(0); 162 } 163 164 /*@C 165 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 166 167 Not Collective 168 169 Input Parameter: 170 . da - the DMDA object 171 172 Output Parameters: 173 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 174 175 Level: intermediate 176 177 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 178 @*/ 179 #undef __FUNCT__ 180 #define __FUNCT__ "DMDAGetElementType" 181 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 182 { 183 DM_DA *dd = (DM_DA*)da->data; 184 PetscFunctionBegin; 185 PetscValidHeaderSpecific(da,DM_CLASSID,1); 186 PetscValidPointer(etype,2); 187 *etype = dd->elementtype; 188 PetscFunctionReturn(0); 189 } 190 191 /*@C 192 DMDAGetElements - Gets an array containing the indices (in local coordinates) 193 of all the local elements 194 195 Not Collective 196 197 Input Parameter: 198 . dm - the DM object 199 200 Output Parameters: 201 + nel - number of local elements 202 . nen - number of element nodes 203 - e - the local indices of the elements' vertices 204 205 Level: intermediate 206 207 .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 208 @*/ 209 #undef __FUNCT__ 210 #define __FUNCT__ "DMDAGetElements" 211 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 212 { 213 DM_DA *da = (DM_DA*)dm->data; 214 PetscErrorCode ierr; 215 PetscFunctionBegin; 216 if (da->dim==-1) { 217 *nel = 0; *nen = 0; *e = PETSC_NULL; 218 } else if (da->dim==1) { 219 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 220 } else if (da->dim==2) { 221 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 222 } else if (da->dim==3) { 223 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 224 } else { 225 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim); 226 } 227 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "DMDARestoreElements" 233 /*@C 234 DMDARestoreElements - Returns an array containing the indices (in local coordinates) 235 of all the local elements obtained with DMDAGetElements() 236 237 Not Collective 238 239 Input Parameter: 240 + dm - the DM object 241 . nel - number of local elements 242 . nen - number of element nodes 243 - e - the local indices of the elements' vertices 244 245 Level: intermediate 246 247 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 248 @*/ 249 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 250 { 251 PetscFunctionBegin; 252 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 253 PetscValidIntPointer(nel,2); 254 PetscValidIntPointer(nen,3); 255 PetscValidPointer(e,4); 256 /* XXX */ 257 PetscFunctionReturn(0); 258 } 259