1 2 #include <petsc-private/daimpl.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 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 16 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 17 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 18 da->ne = 1*(xe - xs - 1); 19 ierr = PetscMalloc((1 + 2*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 20 for (i=xs; i<xe-1; i++) { 21 da->e[cnt++] = (i-Xs); 22 da->e[cnt++] = (i-Xs+1); 23 } 24 } 25 *nel = da->ne; 26 *nen = 2; 27 *e = da->e; 28 PetscFunctionReturn(0); 29 } 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "DMDAGetElements_2D" 33 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 34 { 35 PetscErrorCode ierr; 36 DM_DA *da = (DM_DA*)dm->data; 37 PetscInt i,xs,xe,Xs,Xe; 38 PetscInt j,ys,ye,Ys,Ye; 39 PetscInt cnt=0, cell[4], ns=2, nn=3; 40 PetscInt c, split[] = {0,1,3, 41 2,3,1}; 42 43 PetscFunctionBegin; 44 if (!da->e) { 45 if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;} 46 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;} 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 = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&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 #undef __FUNCT__ 75 #define __FUNCT__ "DMDAGetElements_3D" 76 static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 77 { 78 PetscErrorCode ierr; 79 DM_DA *da = (DM_DA*)dm->data; 80 PetscInt i,xs,xe,Xs,Xe; 81 PetscInt j,ys,ye,Ys,Ye; 82 PetscInt k,zs,ze,Zs,Ze; 83 PetscInt cnt=0, cell[8], ns=6, nn=4; 84 PetscInt c, split[] = {0,1,3,7, 85 0,1,7,4, 86 1,2,3,7, 87 1,2,7,6, 88 1,4,5,7, 89 1,5,6,7}; 90 91 PetscFunctionBegin; 92 if (!da->e) { 93 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;} 94 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;} 95 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 96 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 97 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 98 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 99 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 100 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 101 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 102 for (k=zs; k<ze-1; k++) { 103 for (j=ys; j<ye-1; j++) { 104 for (i=xs; i<xe-1; i++) { 105 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 106 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 107 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 108 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 109 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 110 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 111 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 112 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 113 if (da->elementtype == DMDA_ELEMENT_P1) { 114 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 115 } 116 if (da->elementtype == DMDA_ELEMENT_Q1) { 117 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 118 } 119 } 120 } 121 } 122 } 123 *nel = da->ne; 124 *nen = nn; 125 *e = da->e; 126 PetscFunctionReturn(0); 127 } 128 129 /*@C 130 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 131 132 Not Collective 133 134 Input Parameter: 135 . da - the DMDA object 136 137 Output Parameters: 138 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 139 140 Level: intermediate 141 142 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 143 @*/ 144 #undef __FUNCT__ 145 #define __FUNCT__ "DMDASetElementType" 146 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 147 { 148 DM_DA *dd = (DM_DA*)da->data; 149 PetscErrorCode ierr; 150 151 PetscFunctionBegin; 152 PetscValidHeaderSpecific(da,DM_CLASSID,1); 153 PetscValidLogicalCollectiveEnum(da,etype,2); 154 if (dd->elementtype != etype) { 155 ierr = PetscFree(dd->e);CHKERRQ(ierr); 156 157 dd->elementtype = etype; 158 dd->ne = 0; 159 dd->e = PETSC_NULL; 160 } 161 PetscFunctionReturn(0); 162 } 163 164 /*@C 165 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 166 167 Not Collective 168 169 Input Parameter: 170 . da - the DMDA object 171 172 Output Parameters: 173 . etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1 174 175 Level: intermediate 176 177 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 178 @*/ 179 #undef __FUNCT__ 180 #define __FUNCT__ "DMDAGetElementType" 181 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 182 { 183 DM_DA *dd = (DM_DA*)da->data; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(da,DM_CLASSID,1); 187 PetscValidPointer(etype,2); 188 *etype = dd->elementtype; 189 PetscFunctionReturn(0); 190 } 191 192 /*@C 193 DMDAGetElements - Gets an array containing the indices (in local coordinates) 194 of all the local elements 195 196 Not Collective 197 198 Input Parameter: 199 . dm - the DM object 200 201 Output Parameters: 202 + nel - number of local elements 203 . nen - number of element nodes 204 - e - the local indices of the elements' vertices 205 206 Level: intermediate 207 208 .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 209 @*/ 210 #undef __FUNCT__ 211 #define __FUNCT__ "DMDAGetElements" 212 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 213 { 214 DM_DA *da = (DM_DA*)dm->data; 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 if (da->dim==-1) { 219 *nel = 0; *nen = 0; *e = PETSC_NULL; 220 } else if (da->dim==1) { 221 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 222 } else if (da->dim==2) { 223 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 224 } else if (da->dim==3) { 225 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 226 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim); 227 PetscFunctionReturn(0); 228 } 229 230 #undef __FUNCT__ 231 #define __FUNCT__ "DMDARestoreElements" 232 /*@C 233 DMDARestoreElements - Returns an array containing the indices (in local coordinates) 234 of all the local elements obtained with DMDAGetElements() 235 236 Not Collective 237 238 Input Parameter: 239 + dm - the DM object 240 . nel - number of local elements 241 . nen - number of element nodes 242 - e - the local indices of the elements' vertices 243 244 Level: intermediate 245 246 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 247 @*/ 248 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 249 { 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 252 PetscValidIntPointer(nel,2); 253 PetscValidIntPointer(nen,3); 254 PetscValidPointer(e,4); 255 /* XXX */ 256 PetscFunctionReturn(0); 257 } 258