Actual source code: dageometry.c
1: #include <petscsf.h>
2: #include <petsc/private/dmdaimpl.h>
4: /*@
5: DMDAConvertToCell - Convert a (i,j,k) location in a `DMDA` to its local cell or vertex number
7: Not Collective
9: Input Parameters:
10: + dm - the `DMDA`
11: - s - a `MatStencil` that provides (i,j,k)
13: Output Parameter:
14: . cell - the local cell or vertext number
16: Level: developer
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()`
21: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetGhostCorners()`
22: @*/
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;
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: }
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;
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;
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));
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: }
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;
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));
149: mxlocal = xe - xs - 1;
150: mylocal = ye - ys - 1;
152: dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal);
153: dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal);
155: PetscCall(VecGetLocalSize(pos, &n));
156: PetscCall(VecGetBlockSize(pos, &bs));
157: npoints = n / bs;
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];
165: coor_p[0] = PetscRealPart(coor[2 * p]);
166: coor_p[1] = PetscRealPart(coor[2 * p + 1]);
168: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
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;
175: for (PetscInt d = 0; d < 2; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]);
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;
182: if (mi[0] == (xe - 1)) mi[0]--;
183: if (mi[1] == (ye - 1)) mi[1]--;
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: }
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;
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));
214: mxlocal = xe - xs - 1;
215: mylocal = ye - ys - 1;
216: mzlocal = ze - zs - 1;
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);
222: PetscCall(VecGetLocalSize(pos, &n));
223: PetscCall(VecGetBlockSize(pos, &bs));
224: npoints = n / bs;
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];
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]);
236: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
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;
245: for (PetscInt d = 0; d < 3; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]);
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;
255: if (mi[0] == (xe - 1)) mi[0]--;
256: if (mi[1] == (ye - 1)) mi[1]--;
257: if (mi[2] == (ze - 1)) mi[2]--;
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: }
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;
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: }
288: PetscCall(VecGetLocalSize(pos, &npoints));
289: PetscCall(VecGetBlockSize(pos, &bs));
290: npoints = npoints / bs;
292: PetscCall(PetscMalloc1(npoints, &cells));
293: PetscCall(ISGetIndices(iscell, &boxCells));
295: for (p = 0; p < npoints; p++) {
296: cells[p].rank = 0;
297: cells[p].index = boxCells[p];
298: }
299: PetscCall(ISRestoreIndices(iscell, &boxCells));
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: }