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 (etype == DMDA_ELEMENT_Q1 && dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Q1 elements require you use a stencil type of DMDA_STENCIL_BOX"); 158 if (dd->elementtype != etype) { 159 ierr = PetscFree(dd->e);CHKERRQ(ierr); 160 161 dd->elementtype = etype; 162 dd->ne = 0; 163 dd->e = NULL; 164 } 165 PetscFunctionReturn(0); 166 } 167 168 /*@C 169 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 170 171 Not Collective 172 173 Input Parameter: 174 . da - the DMDA object 175 176 Output Parameters: 177 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 178 179 Level: intermediate 180 181 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 182 @*/ 183 #undef __FUNCT__ 184 #define __FUNCT__ "DMDAGetElementType" 185 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 186 { 187 DM_DA *dd = (DM_DA*)da->data; 188 189 PetscFunctionBegin; 190 PetscValidHeaderSpecific(da,DM_CLASSID,1); 191 PetscValidPointer(etype,2); 192 *etype = dd->elementtype; 193 PetscFunctionReturn(0); 194 } 195 196 /*@C 197 DMDAGetElements - Gets an array containing the indices (in local coordinates) 198 of all the local elements 199 200 Not Collective 201 202 Input Parameter: 203 . dm - the DM object 204 205 Output Parameters: 206 + nel - number of local elements 207 . nen - number of element nodes 208 - e - the local indices of the elements' vertices 209 210 Level: intermediate 211 212 Notes: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 213 214 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 215 216 .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 217 @*/ 218 #undef __FUNCT__ 219 #define __FUNCT__ "DMDAGetElements" 220 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 221 { 222 PetscInt dim; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 227 if (dim==-1) { 228 *nel = 0; *nen = 0; *e = NULL; 229 } else if (dim==1) { 230 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 231 } else if (dim==2) { 232 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 233 } else if (dim==3) { 234 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 235 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 236 PetscFunctionReturn(0); 237 } 238 239 #undef __FUNCT__ 240 #define __FUNCT__ "DMDARestoreElements" 241 /*@C 242 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 243 244 Not Collective 245 246 Input Parameter: 247 + dm - the DM object 248 . nel - number of local elements 249 . nen - number of element nodes 250 - e - the local indices of the elements' vertices 251 252 Level: intermediate 253 254 Note: You should not access these values after you have called this routine. 255 256 This restore signals the DMDA object that you no longer need access to the array information. 257 258 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 259 @*/ 260 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 261 { 262 PetscFunctionBegin; 263 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 264 PetscValidIntPointer(nel,2); 265 PetscValidIntPointer(nen,3); 266 PetscValidPointer(e,4); 267 *nel = 0; 268 *nen = -1; 269 *e = NULL; 270 PetscFunctionReturn(0); 271 } 272