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