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