1 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 5 { 6 PetscErrorCode ierr; 7 DM_DA *da = (DM_DA*)dm->data; 8 PetscInt i,xs,xe,Xs,Xe; 9 PetscInt cnt=0; 10 11 PetscFunctionBegin; 12 if (!da->e) { 13 PetscInt corners[2]; 14 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 corners[0] = (xs -Xs); 26 corners[1] = (xe-1-Xs); 27 ierr = ISCreateGeneral(PETSC_COMM_SELF,2,corners,PETSC_COPY_VALUES,&da->ecorners);CHKERRQ(ierr); 28 } 29 *nel = da->ne; 30 *nen = 2; 31 *e = da->e; 32 PetscFunctionReturn(0); 33 } 34 35 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 36 { 37 PetscErrorCode ierr; 38 DM_DA *da = (DM_DA*)dm->data; 39 PetscInt i,xs,xe,Xs,Xe; 40 PetscInt j,ys,ye,Ys,Ye; 41 PetscInt cnt=0, cell[4], ns=2, nn=3; 42 PetscInt c, split[] = {0,1,3, 43 2,3,1}; 44 45 PetscFunctionBegin; 46 if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;} 47 if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;} 48 if (!da->e) { 49 PetscInt corners[4]; 50 51 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 52 if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;} 53 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;} 54 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 55 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 56 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 57 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 58 da->ne = ns*(xe - xs - 1)*(ye - ys - 1); 59 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 60 for (j=ys; j<ye-1; j++) { 61 for (i=xs; i<xe-1; i++) { 62 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs); 63 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs); 64 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs); 65 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs); 66 if (da->elementtype == DMDA_ELEMENT_P1) { 67 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 68 } 69 if (da->elementtype == DMDA_ELEMENT_Q1) { 70 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 71 } 72 } 73 } 74 corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs); 75 corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs); 76 corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs); 77 corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs); 78 ierr = ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);CHKERRQ(ierr); 79 } 80 *nel = da->ne; 81 *nen = nn; 82 *e = da->e; 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 87 { 88 PetscErrorCode ierr; 89 DM_DA *da = (DM_DA*)dm->data; 90 PetscInt i,xs,xe,Xs,Xe; 91 PetscInt j,ys,ye,Ys,Ye; 92 PetscInt k,zs,ze,Zs,Ze; 93 PetscInt cnt=0, cell[8], ns=6, nn=4; 94 PetscInt c, split[] = {0,1,3,7, 95 0,1,7,4, 96 1,2,3,7, 97 1,2,7,6, 98 1,4,5,7, 99 1,5,6,7}; 100 101 PetscFunctionBegin; 102 if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;} 103 if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;} 104 if (!da->e) { 105 PetscInt corners[8]; 106 107 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 108 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;} 109 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;} 110 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 111 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 112 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 113 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 114 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 115 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 116 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 117 for (k=zs; k<ze-1; k++) { 118 for (j=ys; j<ye-1; j++) { 119 for (i=xs; i<xe-1; i++) { 120 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 121 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 122 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 123 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 124 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 125 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 126 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 127 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 128 if (da->elementtype == DMDA_ELEMENT_P1) { 129 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 130 } 131 if (da->elementtype == DMDA_ELEMENT_Q1) { 132 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 133 } 134 } 135 } 136 } 137 corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 138 corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 139 corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 140 corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 141 corners[4] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 142 corners[5] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 143 corners[6] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 144 corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 145 ierr = ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);CHKERRQ(ierr); 146 } 147 *nel = da->ne; 148 *nen = nn; 149 *e = da->e; 150 PetscFunctionReturn(0); 151 } 152 153 /*@ 154 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 155 156 Not Collective 157 158 Input Parameter: 159 . da - the DMDA object 160 161 Output Parameters: 162 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 163 164 Level: intermediate 165 166 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 167 @*/ 168 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 169 { 170 DM_DA *dd = (DM_DA*)da->data; 171 PetscErrorCode ierr; 172 PetscBool isda; 173 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(da,DM_CLASSID,1); 176 PetscValidLogicalCollectiveEnum(da,etype,2); 177 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 178 if (!isda) PetscFunctionReturn(0); 179 if (dd->elementtype != etype) { 180 ierr = PetscFree(dd->e);CHKERRQ(ierr); 181 ierr = ISDestroy(&dd->ecorners);CHKERRQ(ierr); 182 183 dd->elementtype = etype; 184 dd->ne = 0; 185 dd->e = NULL; 186 } 187 PetscFunctionReturn(0); 188 } 189 190 /*@ 191 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 192 193 Not Collective 194 195 Input Parameter: 196 . da - the DMDA object 197 198 Output Parameters: 199 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 200 201 Level: intermediate 202 203 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 204 @*/ 205 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 206 { 207 DM_DA *dd = (DM_DA*)da->data; 208 PetscErrorCode ierr; 209 PetscBool isda; 210 211 PetscFunctionBegin; 212 PetscValidHeaderSpecific(da,DM_CLASSID,1); 213 PetscValidPointer(etype,2); 214 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 215 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name); 216 *etype = dd->elementtype; 217 PetscFunctionReturn(0); 218 } 219 220 /*@C 221 DMDAGetElements - Gets an array containing the indices (in local coordinates) 222 of all the local elements 223 224 Not Collective 225 226 Input Parameter: 227 . dm - the DM object 228 229 Output Parameters: 230 + nel - number of local elements 231 . nen - number of element nodes 232 - e - the local indices of the elements' vertices 233 234 Level: intermediate 235 236 Notes: 237 Call DMDARestoreElements() once you have finished accessing the elements. 238 239 Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 240 241 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 242 243 Not supported in Fortran 244 245 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 246 @*/ 247 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 248 { 249 PetscInt dim; 250 PetscErrorCode ierr; 251 DM_DA *dd = (DM_DA*)dm->data; 252 PetscBool isda; 253 254 PetscFunctionBegin; 255 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 256 PetscValidIntPointer(nel,2); 257 PetscValidIntPointer(nen,3); 258 PetscValidPointer(e,4); 259 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);CHKERRQ(ierr); 260 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name); 261 if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElements() requires you use a stencil type of DMDA_STENCIL_BOX"); 262 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 263 if (dim==-1) { 264 *nel = 0; *nen = 0; *e = NULL; 265 } else if (dim==1) { 266 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 267 } else if (dim==2) { 268 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 269 } else if (dim==3) { 270 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 271 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 272 PetscFunctionReturn(0); 273 } 274 275 /*@ 276 DMDAGetElementsCornersIS - Gets an index set containing the corner indices (in local coordinates) 277 of the non-overlapping decomposition identified by DMDAGetElements 278 279 Not Collective 280 281 Input Parameter: 282 . dm - the DM object 283 284 Output Parameters: 285 . is - the index set 286 287 Level: intermediate 288 289 Notes: Call DMDARestoreElementsCornersIS() once you have finished accessing the index set. 290 291 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS() 292 @*/ 293 PetscErrorCode DMDAGetElementsCornersIS(DM dm,IS *is) 294 { 295 PetscErrorCode ierr; 296 DM_DA *dd = (DM_DA*)dm->data; 297 PetscBool isda; 298 299 PetscFunctionBegin; 300 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 301 PetscValidPointer(is,2); 302 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);CHKERRQ(ierr); 303 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name); 304 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"); 305 if (!dd->ecorners) { /* compute elements if not yet done */ 306 const PetscInt *e; 307 PetscInt nel,nen; 308 309 ierr = DMDAGetElements(dm,&nel,&nen,&e);CHKERRQ(ierr); 310 ierr = DMDARestoreElements(dm,&nel,&nen,&e);CHKERRQ(ierr); 311 } 312 *is = dd->ecorners; 313 PetscFunctionReturn(0); 314 } 315 316 /*@C 317 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 318 319 Not Collective 320 321 Input Parameter: 322 + dm - the DM object 323 . nel - number of local elements 324 . nen - number of element nodes 325 - e - the local indices of the elements' vertices 326 327 Level: intermediate 328 329 Note: You should not access these values after you have called this routine. 330 331 This restore signals the DMDA object that you no longer need access to the array information. 332 333 Not supported in Fortran 334 335 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 336 @*/ 337 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 338 { 339 PetscFunctionBegin; 340 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 341 PetscValidIntPointer(nel,2); 342 PetscValidIntPointer(nen,3); 343 PetscValidPointer(e,4); 344 *nel = 0; 345 *nen = -1; 346 *e = NULL; 347 PetscFunctionReturn(0); 348 } 349 350 /*@ 351 DMDARestoreElementsCornersIS - Restores the IS obtained with DMDAGetElementsCornersIS() 352 353 Not Collective 354 355 Input Parameter: 356 + dm - the DM object 357 - is - the index set 358 359 Level: intermediate 360 361 Note: 362 363 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElementsCornersIS() 364 @*/ 365 PetscErrorCode DMDARestoreElementsCornersIS(DM dm,IS *is) 366 { 367 PetscFunctionBegin; 368 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 369 PetscValidHeaderSpecific(*is,IS_CLASSID,2); 370 *is = NULL; 371 PetscFunctionReturn(0); 372 } 373