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