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 /*@C 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 /*@C 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 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 213 @*/ 214 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 215 { 216 PetscInt dim; 217 PetscErrorCode ierr; 218 DM_DA *dd = (DM_DA*)dm->data; 219 220 PetscFunctionBegin; 221 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"); 222 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 223 if (dim==-1) { 224 *nel = 0; *nen = 0; *e = NULL; 225 } else if (dim==1) { 226 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 227 } else if (dim==2) { 228 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 229 } else if (dim==3) { 230 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 231 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 232 PetscFunctionReturn(0); 233 } 234 235 /*@C 236 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 237 238 Not Collective 239 240 Input Parameter: 241 + dm - the DM object 242 . nel - number of local elements 243 . nen - number of element nodes 244 - e - the local indices of the elements' vertices 245 246 Level: intermediate 247 248 Note: You should not access these values after you have called this routine. 249 250 This restore signals the DMDA object that you no longer need access to the array information. 251 252 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 253 @*/ 254 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 255 { 256 PetscFunctionBegin; 257 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 258 PetscValidIntPointer(nel,2); 259 PetscValidIntPointer(nen,3); 260 PetscValidPointer(e,4); 261 *nel = 0; 262 *nen = -1; 263 *e = NULL; 264 PetscFunctionReturn(0); 265 } 266