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 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 14 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 15 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 16 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 17 da->ne = 1*(xe - xs - 1); 18 ierr = PetscMalloc1(1 + 2*da->ne,&da->e);CHKERRQ(ierr); 19 for (i=xs; i<xe-1; i++) { 20 da->e[cnt++] = (i-Xs); 21 da->e[cnt++] = (i-Xs+1); 22 } 23 } 24 *nel = da->ne; 25 *nen = 2; 26 *e = da->e; 27 PetscFunctionReturn(0); 28 } 29 30 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 31 { 32 PetscErrorCode ierr; 33 DM_DA *da = (DM_DA*)dm->data; 34 PetscInt i,xs,xe,Xs,Xe; 35 PetscInt j,ys,ye,Ys,Ye; 36 PetscInt cnt=0, cell[4], ns=2, nn=3; 37 PetscInt c, split[] = {0,1,3, 38 2,3,1}; 39 40 PetscFunctionBegin; 41 if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;} 42 if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;} 43 if (!da->e) { 44 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 45 if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;} 46 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;} 47 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 48 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 49 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 50 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 51 da->ne = ns*(xe - xs - 1)*(ye - ys - 1); 52 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 53 for (j=ys; j<ye-1; j++) { 54 for (i=xs; i<xe-1; i++) { 55 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs); 56 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs); 57 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs); 58 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs); 59 if (da->elementtype == DMDA_ELEMENT_P1) { 60 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 61 } 62 if (da->elementtype == DMDA_ELEMENT_Q1) { 63 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 64 } 65 } 66 } 67 } 68 *nel = da->ne; 69 *nen = nn; 70 *e = da->e; 71 PetscFunctionReturn(0); 72 } 73 74 static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 75 { 76 PetscErrorCode ierr; 77 DM_DA *da = (DM_DA*)dm->data; 78 PetscInt i,xs,xe,Xs,Xe; 79 PetscInt j,ys,ye,Ys,Ye; 80 PetscInt k,zs,ze,Zs,Ze; 81 PetscInt cnt=0, cell[8], ns=6, nn=4; 82 PetscInt c, split[] = {0,1,3,7, 83 0,1,7,4, 84 1,2,3,7, 85 1,2,7,6, 86 1,4,5,7, 87 1,5,6,7}; 88 89 PetscFunctionBegin; 90 if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;} 91 if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;} 92 if (!da->e) { 93 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 94 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;} 95 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;} 96 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 97 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 98 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 99 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 100 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 101 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 102 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 103 for (k=zs; k<ze-1; k++) { 104 for (j=ys; j<ye-1; j++) { 105 for (i=xs; i<xe-1; i++) { 106 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 107 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 108 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 109 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 110 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 111 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 112 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 113 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 114 if (da->elementtype == DMDA_ELEMENT_P1) { 115 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 116 } 117 if (da->elementtype == DMDA_ELEMENT_Q1) { 118 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 119 } 120 } 121 } 122 } 123 } 124 *nel = da->ne; 125 *nen = nn; 126 *e = da->e; 127 PetscFunctionReturn(0); 128 } 129 130 /*@ 131 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 132 133 Not Collective 134 135 Input Parameter: 136 . da - the DMDA object 137 138 Output Parameters: 139 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 140 141 Level: intermediate 142 143 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 144 @*/ 145 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 146 { 147 DM_DA *dd = (DM_DA*)da->data; 148 PetscErrorCode ierr; 149 150 PetscFunctionBegin; 151 PetscValidHeaderSpecific(da,DM_CLASSID,1); 152 PetscValidLogicalCollectiveEnum(da,etype,2); 153 if (dd->elementtype != etype) { 154 ierr = PetscFree(dd->e);CHKERRQ(ierr); 155 156 dd->elementtype = etype; 157 dd->ne = 0; 158 dd->e = NULL; 159 } 160 PetscFunctionReturn(0); 161 } 162 163 /*@ 164 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 165 166 Not Collective 167 168 Input Parameter: 169 . da - the DMDA object 170 171 Output Parameters: 172 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 173 174 Level: intermediate 175 176 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 177 @*/ 178 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 179 { 180 DM_DA *dd = (DM_DA*)da->data; 181 182 PetscFunctionBegin; 183 PetscValidHeaderSpecific(da,DM_CLASSID,1); 184 PetscValidPointer(etype,2); 185 *etype = dd->elementtype; 186 PetscFunctionReturn(0); 187 } 188 189 /*@C 190 DMDAGetElements - Gets an array containing the indices (in local coordinates) 191 of all the local elements 192 193 Not Collective 194 195 Input Parameter: 196 . dm - the DM object 197 198 Output Parameters: 199 + nel - number of local elements 200 . nen - number of element nodes 201 - e - the local indices of the elements' vertices 202 203 Level: intermediate 204 205 Notes: 206 Call DMDARestoreElements() once you have finished accessing the elements. 207 208 Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 209 210 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 211 212 Not supported in Fortran 213 214 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 215 @*/ 216 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 217 { 218 PetscInt dim; 219 PetscErrorCode ierr; 220 DM_DA *dd = (DM_DA*)dm->data; 221 222 PetscFunctionBegin; 223 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"); 224 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 225 if (dim==-1) { 226 *nel = 0; *nen = 0; *e = NULL; 227 } else if (dim==1) { 228 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 229 } else if (dim==2) { 230 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 231 } else if (dim==3) { 232 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 233 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 234 PetscFunctionReturn(0); 235 } 236 237 /*@C 238 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 239 240 Not Collective 241 242 Input Parameter: 243 + dm - the DM object 244 . nel - number of local elements 245 . nen - number of element nodes 246 - e - the local indices of the elements' vertices 247 248 Level: intermediate 249 250 Note: You should not access these values after you have called this routine. 251 252 This restore signals the DMDA object that you no longer need access to the array information. 253 254 Not supported in Fortran 255 256 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 257 @*/ 258 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 259 { 260 PetscFunctionBegin; 261 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 262 PetscValidIntPointer(nel,2); 263 PetscValidIntPointer(nen,3); 264 PetscValidPointer(e,4); 265 *nel = 0; 266 *nen = -1; 267 *e = NULL; 268 PetscFunctionReturn(0); 269 } 270