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; 50 Xe += Xs; 51 ye += ys; 52 Ye += Ys; 53 if (xs != Xs && Xs >= 0) xs -= 1; 54 if (ys != Ys && Ys >= 0) ys -= 1; 55 56 PetscCall(DMGetCoordinatesLocal(dmregular, &coor)); 57 PetscCall(VecGetArrayRead(coor, &_coor)); 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 PetscCall(VecRestoreArrayRead(coor, &_coor)); 67 68 mxlocal = xe - xs - 1; 69 mylocal = ye - ys - 1; 70 71 dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal); 72 dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal); 73 74 PetscCall(DMGetBoundingBox(dmregular, gmin, gmax)); 75 76 PetscCall(VecGetLocalSize(pos, &n)); 77 PetscCall(VecGetBlockSize(pos, &bs)); 78 npoints = n / bs; 79 80 PetscCall(PetscMalloc1(npoints, &cellidx)); 81 PetscCall(VecGetArrayRead(pos, &_coor)); 82 for (p = 0; p < npoints; p++) { 83 PetscReal coor_p[2]; 84 PetscInt mi[2]; 85 86 coor_p[0] = PetscRealPart(_coor[2 * p]); 87 coor_p[1] = PetscRealPart(_coor[2 * p + 1]); 88 89 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 90 91 if (coor_p[0] < gmin_l[0]) continue; 92 if (coor_p[0] > gmax_l[0]) continue; 93 if (coor_p[1] < gmin_l[1]) continue; 94 if (coor_p[1] > gmax_l[1]) continue; 95 96 for (d = 0; d < 2; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]); 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 PetscCall(VecRestoreArrayRead(pos, &_coor)); 109 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell)); 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 126 PetscFunctionBegin; 127 PetscCall(DMDAGetCorners(dmregular, &xs, &ys, &zs, &xe, &ye, &ze)); 128 PetscCall(DMDAGetGhostCorners(dmregular, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze)); 129 xe += xs; 130 Xe += Xs; 131 ye += ys; 132 Ye += Ys; 133 ze += zs; 134 Ze += Zs; 135 if (xs != Xs && Xs >= 0) xs -= 1; 136 if (ys != Ys && Ys >= 0) ys -= 1; 137 if (zs != Zs && Zs >= 0) zs -= 1; 138 139 PetscCall(DMGetCoordinatesLocal(dmregular, &coor)); 140 PetscCall(VecGetArrayRead(coor, &_coor)); 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 PetscCall(VecRestoreArrayRead(coor, &_coor)); 152 153 if (xs != Xs) xs -= 1; 154 if (ys != Ys) ys -= 1; 155 if (zs != Zs) zs -= 1; 156 157 mxlocal = xe - xs - 1; 158 mylocal = ye - ys - 1; 159 mzlocal = ze - zs - 1; 160 161 dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal); 162 dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal); 163 dx[2] = (gmax_l[2] - gmin_l[2]) / ((PetscReal)mzlocal); 164 165 PetscCall(DMGetBoundingBox(dmregular, gmin, gmax)); 166 167 PetscCall(VecGetLocalSize(pos, &n)); 168 PetscCall(VecGetBlockSize(pos, &bs)); 169 npoints = n / bs; 170 171 PetscCall(PetscMalloc1(npoints, &cellidx)); 172 PetscCall(VecGetArrayRead(pos, &_coor)); 173 for (p = 0; p < npoints; p++) { 174 PetscReal coor_p[3]; 175 PetscInt mi[3]; 176 177 coor_p[0] = PetscRealPart(_coor[3 * p]); 178 coor_p[1] = PetscRealPart(_coor[3 * p + 1]); 179 coor_p[2] = PetscRealPart(_coor[3 * p + 2]); 180 181 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 182 183 if (coor_p[0] < gmin_l[0]) continue; 184 if (coor_p[0] > gmax_l[0]) continue; 185 if (coor_p[1] < gmin_l[1]) continue; 186 if (coor_p[1] > gmax_l[1]) continue; 187 if (coor_p[2] < gmin_l[2]) continue; 188 if (coor_p[2] > gmax_l[2]) continue; 189 190 for (d = 0; d < 3; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]); 191 192 if (mi[0] < xs) continue; 193 if (mi[0] > (xe - 1)) continue; 194 if (mi[1] < ys) continue; 195 if (mi[1] > (ye - 1)) continue; 196 if (mi[2] < zs) continue; 197 if (mi[2] > (ze - 1)) continue; 198 199 if (mi[0] == (xe - 1)) mi[0]--; 200 if (mi[1] == (ye - 1)) mi[1]--; 201 if (mi[2] == (ze - 1)) mi[2]--; 202 203 cellidx[p] = (mi[0] - xs) + (mi[1] - ys) * mxlocal + (mi[2] - zs) * mxlocal * mylocal; 204 } 205 PetscCall(VecRestoreArrayRead(pos, &_coor)); 206 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell)); 207 PetscFunctionReturn(0); 208 } 209 210 PetscErrorCode DMLocatePoints_DA_Regular(DM dm, Vec pos, DMPointLocationType ltype, PetscSF cellSF) 211 { 212 IS iscell; 213 PetscSFNode *cells; 214 PetscInt p, bs, dim, npoints, nfound; 215 const PetscInt *boxCells; 216 217 PetscFunctionBegin; 218 PetscCall(VecGetBlockSize(pos, &dim)); 219 switch (dim) { 220 case 1: 221 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support not provided for 1D"); 222 case 2: 223 PetscCall(private_DMDALocatePointsIS_2D_Regular(dm, pos, &iscell)); 224 break; 225 case 3: 226 PetscCall(private_DMDALocatePointsIS_3D_Regular(dm, pos, &iscell)); 227 break; 228 default: 229 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupport spatial dimension"); 230 } 231 232 PetscCall(VecGetLocalSize(pos, &npoints)); 233 PetscCall(VecGetBlockSize(pos, &bs)); 234 npoints = npoints / bs; 235 236 PetscCall(PetscMalloc1(npoints, &cells)); 237 PetscCall(ISGetIndices(iscell, &boxCells)); 238 239 for (p = 0; p < npoints; p++) { 240 cells[p].rank = 0; 241 cells[p].index = boxCells[p]; 242 } 243 PetscCall(ISRestoreIndices(iscell, &boxCells)); 244 245 nfound = npoints; 246 PetscCall(PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER)); 247 PetscCall(ISDestroy(&iscell)); 248 PetscFunctionReturn(0); 249 } 250