1 2 #include <petsc-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 13 PetscFunctionBegin; 14 if (!da->e) { 15 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 16 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 17 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 18 da->ne = 1*(xe - xs - 1); 19 ierr = PetscMalloc((1 + 2*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 29 } 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "DMDAGetElements_2D" 33 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 34 { 35 PetscErrorCode ierr; 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}; 42 43 PetscFunctionBegin; 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 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 48 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 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 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 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++) 61 da->e[cnt++] = cell[split[c]]; 62 } 63 if (da->elementtype == DMDA_ELEMENT_Q1) { 64 for (c=0; c<ns*nn; c++) 65 da->e[cnt++] = cell[c]; 66 } 67 } 68 } 69 } 70 *nel = da->ne; 71 *nen = nn; 72 *e = da->e; 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "DMDAGetElements_3D" 78 static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 79 { 80 PetscErrorCode ierr; 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}; 92 93 PetscFunctionBegin; 94 if (!da->e) { 95 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;} 96 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;} 97 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 98 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 99 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 100 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 101 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 102 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 103 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 104 for (k=zs; k<ze-1; k++) { 105 for (j=ys; j<ye-1; j++) { 106 for (i=xs; i<xe-1; i++) { 107 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 108 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 109 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 110 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 111 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 112 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 113 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 114 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 115 if (da->elementtype == DMDA_ELEMENT_P1) { 116 for (c=0; c<ns*nn; c++) 117 da->e[cnt++] = cell[split[c]]; 118 } 119 if (da->elementtype == DMDA_ELEMENT_Q1) { 120 for (c=0; c<ns*nn; c++) 121 da->e[cnt++] = cell[c]; 122 } 123 } 124 } 125 } 126 } 127 *nel = da->ne; 128 *nen = nn; 129 *e = da->e; 130 PetscFunctionReturn(0); 131 } 132 133 /*@C 134 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 135 136 Not Collective 137 138 Input Parameter: 139 . da - the DMDA object 140 141 Output Parameters: 142 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 143 144 Level: intermediate 145 146 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 147 @*/ 148 #undef __FUNCT__ 149 #define __FUNCT__ "DMDASetElementType" 150 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 151 { 152 DM_DA *dd = (DM_DA*)da->data; 153 PetscErrorCode ierr; 154 155 PetscFunctionBegin; 156 PetscValidHeaderSpecific(da,DM_CLASSID,1); 157 PetscValidLogicalCollectiveEnum(da,etype,2); 158 if (dd->elementtype != etype) { 159 ierr = PetscFree(dd->e);CHKERRQ(ierr); 160 dd->elementtype = etype; 161 dd->ne = 0; 162 dd->e = PETSC_NULL; 163 } 164 PetscFunctionReturn(0); 165 } 166 167 /*@C 168 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 169 170 Not Collective 171 172 Input Parameter: 173 . da - the DMDA object 174 175 Output Parameters: 176 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 177 178 Level: intermediate 179 180 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 181 @*/ 182 #undef __FUNCT__ 183 #define __FUNCT__ "DMDAGetElementType" 184 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 185 { 186 DM_DA *dd = (DM_DA*)da->data; 187 188 PetscFunctionBegin; 189 PetscValidHeaderSpecific(da,DM_CLASSID,1); 190 PetscValidPointer(etype,2); 191 *etype = dd->elementtype; 192 PetscFunctionReturn(0); 193 } 194 195 /*@C 196 DMDAGetElements - Gets an array containing the indices (in local coordinates) 197 of all the local elements 198 199 Not Collective 200 201 Input Parameter: 202 . dm - the DM object 203 204 Output Parameters: 205 + nel - number of local elements 206 . nen - number of element nodes 207 - e - the local indices of the elements' vertices 208 209 Level: intermediate 210 211 .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 212 @*/ 213 #undef __FUNCT__ 214 #define __FUNCT__ "DMDAGetElements" 215 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 216 { 217 DM_DA *da = (DM_DA*)dm->data; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 if (da->dim==-1) { 222 *nel = 0; *nen = 0; *e = PETSC_NULL; 223 } else if (da->dim==1) { 224 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 225 } else if (da->dim==2) { 226 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 227 } else if (da->dim==3) { 228 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 229 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim); 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "DMDARestoreElements" 235 /*@C 236 DMDARestoreElements - Returns an array containing the indices (in local coordinates) 237 of all the local elements obtained with DMDAGetElements() 238 239 Not Collective 240 241 Input Parameter: 242 + dm - the DM object 243 . nel - number of local elements 244 . nen - number of element nodes 245 - e - the local indices of the elements' vertices 246 247 Level: intermediate 248 249 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 250 @*/ 251 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 252 { 253 PetscFunctionBegin; 254 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 255 PetscValidIntPointer(nel,2); 256 PetscValidIntPointer(nen,3); 257 PetscValidPointer(e,4); 258 /* XXX */ 259 PetscFunctionReturn(0); 260 } 261