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