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 224 PetscFunctionBegin; 225 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 226 if (dim==-1) { 227 *nel = 0; *nen = 0; *e = NULL; 228 } else if (dim==1) { 229 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 230 } else if (dim==2) { 231 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 232 } else if (dim==3) { 233 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 234 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 235 PetscFunctionReturn(0); 236 } 237 238 #undef __FUNCT__ 239 #define __FUNCT__ "DMDARestoreElements" 240 /*@C 241 DMDARestoreElements - Restores the array 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 Note: You should not access these values after you have called this routine. 254 255 This restore signals the DMDA object that you no longer need access to the array information. 256 257 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 258 @*/ 259 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 260 { 261 PetscFunctionBegin; 262 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 263 PetscValidIntPointer(nel,2); 264 PetscValidIntPointer(nen,3); 265 PetscValidPointer(e,4); 266 *nel = 0; 267 *nen = -1; 268 *e = NULL; 269 PetscFunctionReturn(0); 270 } 271