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