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 distributed array
11: - s - A `MatStencil` that provides (i,j,k)
13: Output Parameter:
14: . cell - the local cell or vertext number
16: Level: developer
18: .seealso: [](sec_struct), `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;
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: }
36: static 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;
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;
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);
63: gmin_l[0] = PetscRealPart(_coor[2 * c0 + 0]);
64: gmin_l[1] = PetscRealPart(_coor[2 * c0 + 1]);
66: gmax_l[0] = PetscRealPart(_coor[2 * c1 + 0]);
67: gmax_l[1] = PetscRealPart(_coor[2 * c1 + 1]);
68: PetscCall(VecRestoreArrayRead(coor, &_coor));
70: mxlocal = xe - xs - 1;
71: mylocal = ye - ys - 1;
73: dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal);
74: dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal);
76: PetscCall(DMGetBoundingBox(dmregular, gmin, gmax));
78: PetscCall(VecGetLocalSize(pos, &n));
79: PetscCall(VecGetBlockSize(pos, &bs));
80: npoints = n / bs;
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];
88: coor_p[0] = PetscRealPart(_coor[2 * p]);
89: coor_p[1] = PetscRealPart(_coor[2 * p + 1]);
91: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
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;
98: for (d = 0; d < 2; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]);
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;
105: if (mi[0] == (xe - 1)) mi[0]--;
106: if (mi[1] == (ye - 1)) mi[1]--;
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: }
115: static 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;
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;
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);
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]);
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));
155: mxlocal = xe - xs - 1;
156: mylocal = ye - ys - 1;
157: mzlocal = ze - zs - 1;
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);
163: PetscCall(DMGetBoundingBox(dmregular, gmin, gmax));
165: PetscCall(VecGetLocalSize(pos, &n));
166: PetscCall(VecGetBlockSize(pos, &bs));
167: npoints = n / bs;
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];
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]);
179: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
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;
188: for (d = 0; d < 3; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]);
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;
197: if (mi[0] == (xe - 1)) mi[0]--;
198: if (mi[1] == (ye - 1)) mi[1]--;
199: if (mi[2] == (ze - 1)) mi[2]--;
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: }
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;
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: }
230: PetscCall(VecGetLocalSize(pos, &npoints));
231: PetscCall(VecGetBlockSize(pos, &bs));
232: npoints = npoints / bs;
234: PetscCall(PetscMalloc1(npoints, &cells));
235: PetscCall(ISGetIndices(iscell, &boxCells));
237: for (p = 0; p < npoints; p++) {
238: cells[p].rank = 0;
239: cells[p].index = boxCells[p];
240: }
241: PetscCall(ISRestoreIndices(iscell, &boxCells));
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: }