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