1 #include <petscsf.h> 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 5 /*@ 6 DMDAConvertToCell - Convert (i,j,k) to local cell number 7 8 Not Collective 9 10 Input Parameter: 11 + da - the distributed array 12 - s - A MatStencil giving (i,j,k) 13 14 Output Parameter: 15 . cell - the local cell number 16 17 Level: developer 18 19 .seealso: DMDAVecGetClosure() 20 @*/ 21 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 22 { 23 DM_DA *da = (DM_DA*) dm->data; 24 const PetscInt dim = dm->dim; 25 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/; 26 const PetscInt il = s.i - da->Xs/da->w, jl = dim > 1 ? s.j - da->Ys : 0, kl = dim > 2 ? s.k - da->Zs : 0; 27 28 PetscFunctionBegin; 29 *cell = -1; 30 if ((s.i < da->Xs/da->w) || (s.i >= da->Xe/da->w)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil i %D should be in [%D, %D)", s.i, da->Xs/da->w, da->Xe/da->w); 31 if ((dim > 1) && ((s.j < da->Ys) || (s.j >= da->Ye))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil j %D should be in [%D, %D)", s.j, da->Ys, da->Ye); 32 if ((dim > 2) && ((s.k < da->Zs) || (s.k >= da->Ze))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil k %D should be in [%D, %D)", s.k, da->Zs, da->Ze); 33 *cell = (kl*my + jl)*mx + il; 34 PetscFunctionReturn(0); 35 } 36 37 PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular,Vec pos,IS *iscell) 38 { 39 PetscInt n,bs,p,npoints; 40 PetscInt xs,xe,Xs,Xe,mxlocal; 41 PetscInt ys,ye,Ys,Ye,mylocal; 42 PetscInt d,c0,c1; 43 PetscReal gmin_l[2],gmax_l[2],dx[2]; 44 PetscReal gmin[2],gmax[2]; 45 PetscInt *cellidx; 46 Vec coor; 47 const PetscScalar *_coor; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 ierr = DMDAGetCorners(dmregular,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 52 ierr = DMDAGetGhostCorners(dmregular,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 53 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 54 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 55 56 mxlocal = xe - xs - 1; 57 mylocal = ye - ys - 1; 58 59 ierr = DMGetCoordinatesLocal(dmregular,&coor);CHKERRQ(ierr); 60 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr); 61 c0 = (xs-Xs) + (ys-Ys)*(Xe-Xs); 62 c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs); 63 64 gmin_l[0] = PetscRealPart(_coor[2*c0+0]); 65 gmin_l[1] = PetscRealPart(_coor[2*c0+1]); 66 67 gmax_l[0] = PetscRealPart(_coor[2*c1+0]); 68 gmax_l[1] = PetscRealPart(_coor[2*c1+1]); 69 70 dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal); 71 dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal); 72 73 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr); 74 75 ierr = DMGetBoundingBox(dmregular,gmin,gmax);CHKERRQ(ierr); 76 77 ierr = VecGetLocalSize(pos,&n);CHKERRQ(ierr); 78 ierr = VecGetBlockSize(pos,&bs);CHKERRQ(ierr); 79 npoints = n/bs; 80 81 ierr = PetscMalloc1(npoints,&cellidx);CHKERRQ(ierr); 82 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 83 for (p=0; p<npoints; p++) { 84 PetscReal coor_p[2]; 85 PetscInt mi[2]; 86 87 coor_p[0] = PetscRealPart(_coor[2*p]); 88 coor_p[1] = PetscRealPart(_coor[2*p+1]); 89 90 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 91 92 if (coor_p[0] < gmin_l[0]) { continue; } 93 if (coor_p[0] > gmax_l[0]) { continue; } 94 if (coor_p[1] < gmin_l[1]) { continue; } 95 if (coor_p[1] > gmax_l[1]) { continue; } 96 97 for (d=0; d<2; d++) { 98 mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] ); 99 } 100 101 if (mi[0] < xs) { continue; } 102 if (mi[0] > (xe-1)) { continue; } 103 if (mi[1] < ys) { continue; } 104 if (mi[1] > (ye-1)) { continue; } 105 106 if (mi[0] == (xe-1)) { mi[0]--; } 107 if (mi[1] == (ye-1)) { mi[1]--; } 108 109 cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal; 110 } 111 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 112 ierr = ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular,Vec pos,IS *iscell) 117 { 118 PetscInt n,bs,p,npoints; 119 PetscInt xs,xe,Xs,Xe,mxlocal; 120 PetscInt ys,ye,Ys,Ye,mylocal; 121 PetscInt zs,ze,Zs,Ze,mzlocal; 122 PetscInt d,c0,c1; 123 PetscReal gmin_l[3],gmax_l[3],dx[3]; 124 PetscReal gmin[3],gmax[3]; 125 PetscInt *cellidx; 126 Vec coor; 127 const PetscScalar *_coor; 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 ierr = DMDAGetCorners(dmregular,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 132 ierr = DMDAGetGhostCorners(dmregular,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 133 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 134 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 135 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 136 137 mxlocal = xe - xs - 1; 138 mylocal = ye - ys - 1; 139 mzlocal = ze - zs - 1; 140 141 ierr = DMGetCoordinatesLocal(dmregular,&coor);CHKERRQ(ierr); 142 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr); 143 c0 = (xs-Xs) + (ys-Ys) *(Xe-Xs) + (zs-Zs) *(Xe-Xs)*(Ye-Ys); 144 c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs) + (ze-2-Zs+1)*(Xe-Xs)*(Ye-Ys); 145 146 gmin_l[0] = PetscRealPart(_coor[3*c0+0]); 147 gmin_l[1] = PetscRealPart(_coor[3*c0+1]); 148 gmin_l[2] = PetscRealPart(_coor[3*c0+2]); 149 150 gmax_l[0] = PetscRealPart(_coor[3*c1+0]); 151 gmax_l[1] = PetscRealPart(_coor[3*c1+1]); 152 gmax_l[2] = PetscRealPart(_coor[3*c1+2]); 153 154 dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal); 155 dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal); 156 dx[2] = (gmax_l[2]-gmin_l[2])/((PetscReal)mzlocal); 157 158 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr); 159 160 ierr = DMGetBoundingBox(dmregular,gmin,gmax);CHKERRQ(ierr); 161 162 ierr = VecGetLocalSize(pos,&n);CHKERRQ(ierr); 163 ierr = VecGetBlockSize(pos,&bs);CHKERRQ(ierr); 164 npoints = n/bs; 165 166 ierr = PetscMalloc1(npoints,&cellidx);CHKERRQ(ierr); 167 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 168 for (p=0; p<npoints; p++) { 169 PetscReal coor_p[3]; 170 PetscInt mi[3]; 171 172 coor_p[0] = PetscRealPart(_coor[3*p]); 173 coor_p[1] = PetscRealPart(_coor[3*p+1]); 174 coor_p[2] = PetscRealPart(_coor[3*p+2]); 175 176 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 177 178 if (coor_p[0] < gmin_l[0]) { continue; } 179 if (coor_p[0] > gmax_l[0]) { continue; } 180 if (coor_p[1] < gmin_l[1]) { continue; } 181 if (coor_p[1] > gmax_l[1]) { continue; } 182 if (coor_p[2] < gmin_l[2]) { continue; } 183 if (coor_p[2] > gmax_l[2]) { continue; } 184 185 for (d=0; d<3; d++) { 186 mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] ); 187 } 188 189 if (mi[0] < xs) { continue; } 190 if (mi[0] > (xe-1)) { continue; } 191 if (mi[1] < ys) { continue; } 192 if (mi[1] > (ye-1)) { continue; } 193 if (mi[2] < zs) { continue; } 194 if (mi[2] > (ze-1)) { continue; } 195 196 if (mi[0] == (xe-1)) { mi[0]--; } 197 if (mi[1] == (ye-1)) { mi[1]--; } 198 if (mi[2] == (ze-1)) { mi[2]--; } 199 200 cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal + (mi[2]-zs) * mxlocal * mylocal; 201 } 202 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 203 ierr = ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 PetscErrorCode DMLocatePoints_DA_Regular(DM dm,Vec pos,DMPointLocationType ltype,PetscSF cellSF) 208 { 209 IS iscell; 210 PetscSFNode *cells; 211 PetscInt p,bs,dim,npoints,nfound; 212 const PetscInt *boxCells; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 ierr = VecGetBlockSize(pos,&dim);CHKERRQ(ierr); 217 switch (dim) { 218 case 1: 219 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support not provided for 1D"); 220 break; 221 case 2: 222 ierr = private_DMDALocatePointsIS_2D_Regular(dm,pos,&iscell);CHKERRQ(ierr); 223 break; 224 case 3: 225 ierr = private_DMDALocatePointsIS_3D_Regular(dm,pos,&iscell);CHKERRQ(ierr); 226 break; 227 default: 228 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupport spatial dimension"); 229 break; 230 } 231 232 ierr = VecGetLocalSize(pos,&npoints);CHKERRQ(ierr); 233 ierr = VecGetBlockSize(pos,&bs);CHKERRQ(ierr); 234 npoints = npoints / bs; 235 236 ierr = PetscMalloc1(npoints, &cells);CHKERRQ(ierr); 237 ierr = ISGetIndices(iscell, &boxCells);CHKERRQ(ierr); 238 239 for (p=0; p<npoints; p++) { 240 cells[p].rank = 0; 241 cells[p].index = boxCells[p]; 242 } 243 ierr = ISRestoreIndices(iscell, &boxCells);CHKERRQ(ierr); 244 245 nfound = npoints; 246 ierr = PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);CHKERRQ(ierr); 247 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250