1 #include <petscsf.h> 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 /*@ 5 DMDAConvertToCell - Convert a (i,j,k) location in a `DMDA` to its local cell or vertex number 6 7 Not Collective 8 9 Input Parameters: 10 + dm - the distributed array 11 - s - A `MatStencil` that provides (i,j,k) 12 13 Output Parameter: 14 . cell - the local cell or vertext number 15 16 Level: developer 17 18 .seealso: [](sec_struct), `DM`, `DMDA` 19 @*/ 20 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 21 { 22 DM_DA *da = (DM_DA *)dm->data; 23 const PetscInt dim = dm->dim; 24 const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/; 25 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; 26 27 PetscFunctionBegin; 28 *cell = -1; 29 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); 30 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); 31 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); 32 *cell = (kl * my + jl) * mx + il; 33 PetscFunctionReturn(PETSC_SUCCESS); 34 } 35 36 PetscErrorCode DMGetLocalBoundingBox_DA(DM da, PetscReal lmin[], PetscReal lmax[], PetscInt cs[], PetscInt ce[]) 37 { 38 PetscInt xs, xe, Xs, Xe; 39 PetscInt ys, ye, Ys, Ye; 40 PetscInt zs, ze, Zs, Ze; 41 PetscInt dim, M, N, P, c0, c1; 42 PetscReal gmax[3] = {0., 0., 0.}; 43 const PetscReal *L, *Lstart; 44 Vec coordinates; 45 const PetscScalar *coor; 46 DMBoundaryType bx, by, bz; 47 48 PetscFunctionBegin; 49 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xe, &ye, &ze)); 50 PetscCall(DMDAGetGhostCorners(da, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze)); 51 PetscCall(DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL)); 52 // Convert from widths to endpoints 53 xe += xs; 54 Xe += Xs; 55 ye += ys; 56 Ye += Ys; 57 ze += zs; 58 Ze += Zs; 59 // What is this doing? 60 if (xs != Xs && Xs >= 0) xs -= 1; 61 if (ys != Ys && Ys >= 0) ys -= 1; 62 if (zs != Zs && Zs >= 0) zs -= 1; 63 64 PetscCall(DMGetCoordinatesLocal(da, &coordinates)); 65 if (!coordinates) { 66 PetscCall(DMGetLocalBoundingIndices_DMDA(da, lmin, lmax)); 67 for (PetscInt d = 0; d < dim; ++d) { 68 if (cs) cs[d] = lmin[d]; 69 if (ce) ce[d] = lmax[d]; 70 } 71 PetscFunctionReturn(PETSC_SUCCESS); 72 } 73 PetscCall(VecGetArrayRead(coordinates, &coor)); 74 switch (dim) { 75 case 1: 76 c0 = (xs - Xs); 77 c1 = (xe - 2 - Xs + 1); 78 break; 79 case 2: 80 c0 = (xs - Xs) + (ys - Ys) * (Xe - Xs); 81 c1 = (xe - 2 - Xs + 1) + (ye - 2 - Ys + 1) * (Xe - Xs); 82 break; 83 case 3: 84 c0 = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys); 85 c1 = (xe - 2 - Xs + 1) + (ye - 2 - Ys + 1) * (Xe - Xs) + (ze - 2 - Zs + 1) * (Xe - Xs) * (Ye - Ys); 86 break; 87 default: 88 SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Invalid dimension %" PetscInt_FMT " for DMDA", dim); 89 } 90 for (PetscInt d = 0; d < dim; ++d) { 91 lmin[d] = PetscRealPart(coor[c0 * dim + d]); 92 lmax[d] = PetscRealPart(coor[c1 * dim + d]); 93 } 94 PetscCall(VecRestoreArrayRead(coordinates, &coor)); 95 96 PetscCall(DMGetPeriodicity(da, NULL, &Lstart, &L)); 97 if (L) { 98 for (PetscInt d = 0; d < dim; ++d) 99 if (L[d] > 0.0) gmax[d] = Lstart[d] + L[d]; 100 } 101 // Must check for periodic boundary 102 if (bx == DM_BOUNDARY_PERIODIC && xe == M) { 103 lmax[0] = gmax[0]; 104 ++xe; 105 } 106 if (by == DM_BOUNDARY_PERIODIC && ye == N) { 107 lmax[1] = gmax[1]; 108 ++ye; 109 } 110 if (bz == DM_BOUNDARY_PERIODIC && ze == P) { 111 lmax[2] = gmax[2]; 112 ++ze; 113 } 114 if (cs) { 115 cs[0] = xs; 116 if (dim > 1) cs[1] = ys; 117 if (dim > 2) cs[2] = zs; 118 } 119 if (ce) { 120 ce[0] = xe; 121 if (dim > 1) ce[1] = ye; 122 if (dim > 2) ce[2] = ze; 123 } 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 static PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular, Vec pos, IS *iscell) 128 { 129 PetscInt n, bs, npoints; 130 PetscInt cs[2], ce[2]; 131 PetscInt xs, xe, mxlocal; 132 PetscInt ys, ye, mylocal; 133 PetscReal gmin_l[2], gmax_l[2], dx[2]; 134 PetscReal gmin[2], gmax[2]; 135 PetscInt *cellidx; 136 const PetscScalar *coor; 137 138 PetscFunctionBegin; 139 PetscCall(DMGetLocalBoundingBox_DA(dmregular, gmin_l, gmax_l, cs, ce)); 140 xs = cs[0]; 141 ys = cs[1]; 142 xe = ce[0]; 143 ye = ce[1]; 144 PetscCall(DMGetBoundingBox(dmregular, gmin, gmax)); 145 146 mxlocal = xe - xs - 1; 147 mylocal = ye - ys - 1; 148 149 dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal); 150 dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal); 151 152 PetscCall(VecGetLocalSize(pos, &n)); 153 PetscCall(VecGetBlockSize(pos, &bs)); 154 npoints = n / bs; 155 156 PetscCall(PetscMalloc1(npoints, &cellidx)); 157 PetscCall(VecGetArrayRead(pos, &coor)); 158 for (PetscInt p = 0; p < npoints; p++) { 159 PetscReal coor_p[2]; 160 PetscInt mi[2]; 161 162 coor_p[0] = PetscRealPart(coor[2 * p]); 163 coor_p[1] = PetscRealPart(coor[2 * p + 1]); 164 165 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 166 167 if (coor_p[0] < gmin_l[0]) continue; 168 if (coor_p[0] > gmax_l[0]) continue; 169 if (coor_p[1] < gmin_l[1]) continue; 170 if (coor_p[1] > gmax_l[1]) continue; 171 172 for (PetscInt d = 0; d < 2; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]); 173 174 if (mi[0] < xs) continue; 175 if (mi[0] > (xe - 1)) continue; 176 if (mi[1] < ys) continue; 177 if (mi[1] > (ye - 1)) continue; 178 179 if (mi[0] == (xe - 1)) mi[0]--; 180 if (mi[1] == (ye - 1)) mi[1]--; 181 182 cellidx[p] = (mi[0] - xs) + (mi[1] - ys) * mxlocal; 183 } 184 PetscCall(VecRestoreArrayRead(pos, &coor)); 185 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell)); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 static PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular, Vec pos, IS *iscell) 190 { 191 PetscInt n, bs, npoints; 192 PetscInt cs[3], ce[3]; 193 PetscInt xs, xe, mxlocal; 194 PetscInt ys, ye, mylocal; 195 PetscInt zs, ze, mzlocal; 196 PetscReal gmin_l[3], gmax_l[3], dx[3]; 197 PetscReal gmin[3], gmax[3]; 198 PetscInt *cellidx; 199 const PetscScalar *coor; 200 201 PetscFunctionBegin; 202 PetscCall(DMGetLocalBoundingBox_DA(dmregular, gmin_l, gmax_l, cs, ce)); 203 xs = cs[0]; 204 ys = cs[1]; 205 zs = cs[2]; 206 xe = ce[0]; 207 ye = ce[1]; 208 ze = ce[2]; 209 PetscCall(DMGetBoundingBox(dmregular, gmin, gmax)); 210 211 mxlocal = xe - xs - 1; 212 mylocal = ye - ys - 1; 213 mzlocal = ze - zs - 1; 214 215 dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal); 216 dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal); 217 dx[2] = (gmax_l[2] - gmin_l[2]) / ((PetscReal)mzlocal); 218 219 PetscCall(VecGetLocalSize(pos, &n)); 220 PetscCall(VecGetBlockSize(pos, &bs)); 221 npoints = n / bs; 222 223 PetscCall(PetscMalloc1(npoints, &cellidx)); 224 PetscCall(VecGetArrayRead(pos, &coor)); 225 for (PetscInt p = 0; p < npoints; p++) { 226 PetscReal coor_p[3]; 227 PetscInt mi[3]; 228 229 coor_p[0] = PetscRealPart(coor[3 * p]); 230 coor_p[1] = PetscRealPart(coor[3 * p + 1]); 231 coor_p[2] = PetscRealPart(coor[3 * p + 2]); 232 233 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 234 235 if (coor_p[0] < gmin_l[0]) continue; 236 if (coor_p[0] > gmax_l[0]) continue; 237 if (coor_p[1] < gmin_l[1]) continue; 238 if (coor_p[1] > gmax_l[1]) continue; 239 if (coor_p[2] < gmin_l[2]) continue; 240 if (coor_p[2] > gmax_l[2]) continue; 241 242 for (PetscInt d = 0; d < 3; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]); 243 244 // TODO: Check for periodicity here 245 if (mi[0] < xs) continue; 246 if (mi[0] > (xe - 1)) continue; 247 if (mi[1] < ys) continue; 248 if (mi[1] > (ye - 1)) continue; 249 if (mi[2] < zs) continue; 250 if (mi[2] > (ze - 1)) continue; 251 252 if (mi[0] == (xe - 1)) mi[0]--; 253 if (mi[1] == (ye - 1)) mi[1]--; 254 if (mi[2] == (ze - 1)) mi[2]--; 255 256 cellidx[p] = (mi[0] - xs) + (mi[1] - ys) * mxlocal + (mi[2] - zs) * mxlocal * mylocal; 257 } 258 PetscCall(VecRestoreArrayRead(pos, &coor)); 259 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell)); 260 PetscFunctionReturn(PETSC_SUCCESS); 261 } 262 263 PetscErrorCode DMLocatePoints_DA_Regular(DM dm, Vec pos, DMPointLocationType ltype, PetscSF cellSF) 264 { 265 IS iscell; 266 PetscSFNode *cells; 267 PetscInt p, bs, dim, npoints, nfound; 268 const PetscInt *boxCells; 269 270 PetscFunctionBegin; 271 PetscCall(VecGetBlockSize(pos, &dim)); 272 switch (dim) { 273 case 1: 274 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support not provided for 1D"); 275 case 2: 276 PetscCall(private_DMDALocatePointsIS_2D_Regular(dm, pos, &iscell)); 277 break; 278 case 3: 279 PetscCall(private_DMDALocatePointsIS_3D_Regular(dm, pos, &iscell)); 280 break; 281 default: 282 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported spatial dimension"); 283 } 284 285 PetscCall(VecGetLocalSize(pos, &npoints)); 286 PetscCall(VecGetBlockSize(pos, &bs)); 287 npoints = npoints / bs; 288 289 PetscCall(PetscMalloc1(npoints, &cells)); 290 PetscCall(ISGetIndices(iscell, &boxCells)); 291 292 for (p = 0; p < npoints; p++) { 293 cells[p].rank = 0; 294 cells[p].index = boxCells[p]; 295 } 296 PetscCall(ISRestoreIndices(iscell, &boxCells)); 297 298 nfound = npoints; 299 PetscCall(PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER)); 300 PetscCall(ISDestroy(&iscell)); 301 PetscFunctionReturn(PETSC_SUCCESS); 302 } 303