1 2 #include <petsc/private/dmdaimpl.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 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 16 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 17 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 18 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 19 da->ne = 1*(xe - xs - 1); 20 ierr = PetscMalloc1(1 + 2*da->ne,&da->e);CHKERRQ(ierr); 21 for (i=xs; i<xe-1; i++) { 22 da->e[cnt++] = (i-Xs); 23 da->e[cnt++] = (i-Xs+1); 24 } 25 } 26 *nel = da->ne; 27 *nen = 2; 28 *e = da->e; 29 PetscFunctionReturn(0); 30 } 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "DMDAGetElements_2D" 34 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 35 { 36 PetscErrorCode ierr; 37 DM_DA *da = (DM_DA*)dm->data; 38 PetscInt i,xs,xe,Xs,Xe; 39 PetscInt j,ys,ye,Ys,Ye; 40 PetscInt cnt=0, cell[4], ns=2, nn=3; 41 PetscInt c, split[] = {0,1,3, 42 2,3,1}; 43 44 PetscFunctionBegin; 45 if (!da->e) { 46 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 47 if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;} 48 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;} 49 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 50 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 51 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 52 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 53 da->ne = ns*(xe - xs - 1)*(ye - ys - 1); 54 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 55 for (j=ys; j<ye-1; j++) { 56 for (i=xs; i<xe-1; i++) { 57 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs); 58 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs); 59 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs); 60 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs); 61 if (da->elementtype == DMDA_ELEMENT_P1) { 62 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 63 } 64 if (da->elementtype == DMDA_ELEMENT_Q1) { 65 for (c=0; c<ns*nn; c++) 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->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 96 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;} 97 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;} 98 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 99 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 100 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 101 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 102 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 103 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 104 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 105 for (k=zs; k<ze-1; k++) { 106 for (j=ys; j<ye-1; j++) { 107 for (i=xs; i<xe-1; i++) { 108 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 109 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 110 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 111 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 112 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 113 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 114 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 115 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 116 if (da->elementtype == DMDA_ELEMENT_P1) { 117 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 118 } 119 if (da->elementtype == DMDA_ELEMENT_Q1) { 120 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 121 } 122 } 123 } 124 } 125 } 126 *nel = da->ne; 127 *nen = nn; 128 *e = da->e; 129 PetscFunctionReturn(0); 130 } 131 132 /*@C 133 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 134 135 Not Collective 136 137 Input Parameter: 138 . da - the DMDA object 139 140 Output Parameters: 141 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 142 143 Level: intermediate 144 145 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 146 @*/ 147 #undef __FUNCT__ 148 #define __FUNCT__ "DMDASetElementType" 149 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 150 { 151 DM_DA *dd = (DM_DA*)da->data; 152 PetscErrorCode ierr; 153 154 PetscFunctionBegin; 155 PetscValidHeaderSpecific(da,DM_CLASSID,1); 156 PetscValidLogicalCollectiveEnum(da,etype,2); 157 if (dd->elementtype != etype) { 158 ierr = PetscFree(dd->e);CHKERRQ(ierr); 159 160 dd->elementtype = etype; 161 dd->ne = 0; 162 dd->e = 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 Notes: 210 Call DMDARestoreElements() once you have finished accessing the elements. 211 212 Level: intermediate 213 214 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 215 @*/ 216 #undef __FUNCT__ 217 #define __FUNCT__ "DMDAGetElements" 218 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 219 { 220 PetscInt dim; 221 PetscErrorCode ierr; 222 223 PetscFunctionBegin; 224 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 225 if (dim==-1) { 226 *nel = 0; *nen = 0; *e = NULL; 227 } else if (dim==1) { 228 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 229 } else if (dim==2) { 230 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 231 } else if (dim==3) { 232 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 233 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 234 PetscFunctionReturn(0); 235 } 236 237 #undef __FUNCT__ 238 #define __FUNCT__ "DMDARestoreElements" 239 /*@C 240 DMDARestoreElements - Restores an array containing the indices (in local coordinates) 241 of all the local elements obtained with DMDAGetElements() 242 243 Not Collective 244 245 Input Parameter: 246 + dm - the DM object 247 . nel - number of local elements 248 . nen - number of element nodes 249 - e - the local indices of the elements' vertices 250 251 Level: intermediate 252 253 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 254 @*/ 255 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 256 { 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 259 PetscValidIntPointer(nel,2); 260 PetscValidIntPointer(nen,3); 261 PetscValidPointer(e,4); 262 /* XXX */ 263 PetscFunctionReturn(0); 264 } 265