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