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 @*/
DMDAConvertToCell(DM dm,MatStencil s,PetscInt * cell)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
DMGetLocalBoundingBox_DA(DM da,PetscReal lmin[],PetscReal lmax[],PetscInt cs[],PetscInt ce[])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
private_DMDALocatePointsIS_2D_Regular(DM dmregular,Vec pos,IS * iscell)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
private_DMDALocatePointsIS_3D_Regular(DM dmregular,Vec pos,IS * iscell)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
DMLocatePoints_DA_Regular(DM dm,Vec pos,DMPointLocationType ltype,PetscSF cellSF)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