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