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 if (xs != Xs) xs -= 1; 52 ye += ys; 53 Ye += Ys; 54 if (ys != Ys) ys -= 1; 55 56 mxlocal = xe - xs - 1; 57 mylocal = ye - ys - 1; 58 59 PetscCall(DMGetCoordinatesLocal(dmregular, &coor)); 60 PetscCall(VecGetArrayRead(coor, &_coor)); 61 c0 = (xs - Xs) + (ys - Ys) * (Xe - Xs); 62 c1 = (xe - 2 - Xs + 1) + (ye - 2 - Ys + 1) * (Xe - Xs); 63 64 gmin_l[0] = PetscRealPart(_coor[2 * c0 + 0]); 65 gmin_l[1] = PetscRealPart(_coor[2 * c0 + 1]); 66 67 gmax_l[0] = PetscRealPart(_coor[2 * c1 + 0]); 68 gmax_l[1] = PetscRealPart(_coor[2 * c1 + 1]); 69 70 dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal); 71 dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal); 72 73 PetscCall(VecRestoreArrayRead(coor, &_coor)); 74 75 PetscCall(DMGetBoundingBox(dmregular, gmin, gmax)); 76 77 PetscCall(VecGetLocalSize(pos, &n)); 78 PetscCall(VecGetBlockSize(pos, &bs)); 79 npoints = n / bs; 80 81 PetscCall(PetscMalloc1(npoints, &cellidx)); 82 PetscCall(VecGetArrayRead(pos, &_coor)); 83 for (p = 0; p < npoints; p++) { 84 PetscReal coor_p[2]; 85 PetscInt mi[2]; 86 87 coor_p[0] = PetscRealPart(_coor[2 * p]); 88 coor_p[1] = PetscRealPart(_coor[2 * p + 1]); 89 90 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 91 92 if (coor_p[0] < gmin_l[0]) continue; 93 if (coor_p[0] > gmax_l[0]) continue; 94 if (coor_p[1] < gmin_l[1]) continue; 95 if (coor_p[1] > gmax_l[1]) continue; 96 97 for (d = 0; d < 2; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]); 98 99 if (mi[0] < xs) continue; 100 if (mi[0] > (xe - 1)) continue; 101 if (mi[1] < ys) continue; 102 if (mi[1] > (ye - 1)) continue; 103 104 if (mi[0] == (xe - 1)) mi[0]--; 105 if (mi[1] == (ye - 1)) mi[1]--; 106 107 cellidx[p] = (mi[0] - xs) + (mi[1] - ys) * mxlocal; 108 } 109 PetscCall(VecRestoreArrayRead(pos, &_coor)); 110 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell)); 111 PetscFunctionReturn(0); 112 } 113 114 PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular, Vec pos, IS *iscell) 115 { 116 PetscInt n, bs, p, npoints; 117 PetscInt xs, xe, Xs, Xe, mxlocal; 118 PetscInt ys, ye, Ys, Ye, mylocal; 119 PetscInt zs, ze, Zs, Ze, mzlocal; 120 PetscInt d, c0, c1; 121 PetscReal gmin_l[3], gmax_l[3], dx[3]; 122 PetscReal gmin[3], gmax[3]; 123 PetscInt *cellidx; 124 Vec coor; 125 const PetscScalar *_coor; 126 127 PetscFunctionBegin; 128 PetscCall(DMDAGetCorners(dmregular, &xs, &ys, &zs, &xe, &ye, &ze)); 129 PetscCall(DMDAGetGhostCorners(dmregular, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze)); 130 xe += xs; 131 Xe += Xs; 132 if (xs != Xs) xs -= 1; 133 ye += ys; 134 Ye += Ys; 135 if (ys != Ys) ys -= 1; 136 ze += zs; 137 Ze += Zs; 138 if (zs != Zs) zs -= 1; 139 140 mxlocal = xe - xs - 1; 141 mylocal = ye - ys - 1; 142 mzlocal = ze - zs - 1; 143 144 PetscCall(DMGetCoordinatesLocal(dmregular, &coor)); 145 PetscCall(VecGetArrayRead(coor, &_coor)); 146 c0 = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys); 147 c1 = (xe - 2 - Xs + 1) + (ye - 2 - Ys + 1) * (Xe - Xs) + (ze - 2 - Zs + 1) * (Xe - Xs) * (Ye - Ys); 148 149 gmin_l[0] = PetscRealPart(_coor[3 * c0 + 0]); 150 gmin_l[1] = PetscRealPart(_coor[3 * c0 + 1]); 151 gmin_l[2] = PetscRealPart(_coor[3 * c0 + 2]); 152 153 gmax_l[0] = PetscRealPart(_coor[3 * c1 + 0]); 154 gmax_l[1] = PetscRealPart(_coor[3 * c1 + 1]); 155 gmax_l[2] = PetscRealPart(_coor[3 * c1 + 2]); 156 157 dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal); 158 dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal); 159 dx[2] = (gmax_l[2] - gmin_l[2]) / ((PetscReal)mzlocal); 160 161 PetscCall(VecRestoreArrayRead(coor, &_coor)); 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(0); 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, "Unsupport 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(0); 247 } 248