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 Level: intermediate 210 211 Notes: 212 Call DMDARestoreElements() once you have finished accessing the elements. 213 214 Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 215 216 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 217 218 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 219 @*/ 220 #undef __FUNCT__ 221 #define __FUNCT__ "DMDAGetElements" 222 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 223 { 224 PetscInt dim; 225 PetscErrorCode ierr; 226 DM_DA *dd = (DM_DA*)dm->data; 227 228 PetscFunctionBegin; 229 if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX"); 230 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 231 if (dim==-1) { 232 *nel = 0; *nen = 0; *e = NULL; 233 } else if (dim==1) { 234 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 235 } else if (dim==2) { 236 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 237 } else if (dim==3) { 238 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 239 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "DMDARestoreElements" 245 /*@C 246 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 247 248 Not Collective 249 250 Input Parameter: 251 + dm - the DM object 252 . nel - number of local elements 253 . nen - number of element nodes 254 - e - the local indices of the elements' vertices 255 256 Level: intermediate 257 258 Note: You should not access these values after you have called this routine. 259 260 This restore signals the DMDA object that you no longer need access to the array information. 261 262 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 263 @*/ 264 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 265 { 266 PetscFunctionBegin; 267 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 268 PetscValidIntPointer(nel,2); 269 PetscValidIntPointer(nen,3); 270 PetscValidPointer(e,4); 271 *nel = 0; 272 *nen = -1; 273 *e = NULL; 274 PetscFunctionReturn(0); 275 } 276