Actual source code: dageometry.c
1: #include <petscsf.h>
2: #include <petsc/private/dmdaimpl.h>
4: /*@
5: DMDAConvertToCell - Convert (i,j,k) to local cell number
7: Not Collective
9: Input Parameters:
10: + da - the distributed array
11: - s - A `MatStencil` giving (i,j,k)
13: Output Parameter:
14: . cell - the local cell number
16: Level: developer
18: .seealso: `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: *cell = -1;
31: *cell = (kl * my + jl) * mx + il;
32: return 0;
33: }
35: PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular, Vec pos, IS *iscell)
36: {
37: PetscInt n, bs, p, npoints;
38: PetscInt xs, xe, Xs, Xe, mxlocal;
39: PetscInt ys, ye, Ys, Ye, mylocal;
40: PetscInt d, c0, c1;
41: PetscReal gmin_l[2], gmax_l[2], dx[2];
42: PetscReal gmin[2], gmax[2];
43: PetscInt *cellidx;
44: Vec coor;
45: const PetscScalar *_coor;
47: DMDAGetCorners(dmregular, &xs, &ys, NULL, &xe, &ye, NULL);
48: DMDAGetGhostCorners(dmregular, &Xs, &Ys, NULL, &Xe, &Ye, NULL);
49: xe += xs;
50: Xe += Xs;
51: ye += ys;
52: Ye += Ys;
53: if (xs != Xs && Xs >= 0) xs -= 1;
54: if (ys != Ys && Ys >= 0) ys -= 1;
56: DMGetCoordinatesLocal(dmregular, &coor);
57: VecGetArrayRead(coor, &_coor);
58: c0 = (xs - Xs) + (ys - Ys) * (Xe - Xs);
59: c1 = (xe - 2 - Xs + 1) + (ye - 2 - Ys + 1) * (Xe - Xs);
61: gmin_l[0] = PetscRealPart(_coor[2 * c0 + 0]);
62: gmin_l[1] = PetscRealPart(_coor[2 * c0 + 1]);
64: gmax_l[0] = PetscRealPart(_coor[2 * c1 + 0]);
65: gmax_l[1] = PetscRealPart(_coor[2 * c1 + 1]);
66: VecRestoreArrayRead(coor, &_coor);
68: mxlocal = xe - xs - 1;
69: mylocal = ye - ys - 1;
71: dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal);
72: dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal);
74: DMGetBoundingBox(dmregular, gmin, gmax);
76: VecGetLocalSize(pos, &n);
77: VecGetBlockSize(pos, &bs);
78: npoints = n / bs;
80: PetscMalloc1(npoints, &cellidx);
81: VecGetArrayRead(pos, &_coor);
82: for (p = 0; p < npoints; p++) {
83: PetscReal coor_p[2];
84: PetscInt mi[2];
86: coor_p[0] = PetscRealPart(_coor[2 * p]);
87: coor_p[1] = PetscRealPart(_coor[2 * p + 1]);
89: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
91: if (coor_p[0] < gmin_l[0]) continue;
92: if (coor_p[0] > gmax_l[0]) continue;
93: if (coor_p[1] < gmin_l[1]) continue;
94: if (coor_p[1] > gmax_l[1]) continue;
96: for (d = 0; d < 2; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]);
98: if (mi[0] < xs) continue;
99: if (mi[0] > (xe - 1)) continue;
100: if (mi[1] < ys) continue;
101: if (mi[1] > (ye - 1)) continue;
103: if (mi[0] == (xe - 1)) mi[0]--;
104: if (mi[1] == (ye - 1)) mi[1]--;
106: cellidx[p] = (mi[0] - xs) + (mi[1] - ys) * mxlocal;
107: }
108: VecRestoreArrayRead(pos, &_coor);
109: ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell);
110: return 0;
111: }
113: PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular, Vec pos, IS *iscell)
114: {
115: PetscInt n, bs, p, npoints;
116: PetscInt xs, xe, Xs, Xe, mxlocal;
117: PetscInt ys, ye, Ys, Ye, mylocal;
118: PetscInt zs, ze, Zs, Ze, mzlocal;
119: PetscInt d, c0, c1;
120: PetscReal gmin_l[3], gmax_l[3], dx[3];
121: PetscReal gmin[3], gmax[3];
122: PetscInt *cellidx;
123: Vec coor;
124: const PetscScalar *_coor;
126: DMDAGetCorners(dmregular, &xs, &ys, &zs, &xe, &ye, &ze);
127: DMDAGetGhostCorners(dmregular, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze);
128: xe += xs;
129: Xe += Xs;
130: ye += ys;
131: Ye += Ys;
132: ze += zs;
133: Ze += Zs;
134: if (xs != Xs && Xs >= 0) xs -= 1;
135: if (ys != Ys && Ys >= 0) ys -= 1;
136: if (zs != Zs && Zs >= 0) zs -= 1;
138: DMGetCoordinatesLocal(dmregular, &coor);
139: VecGetArrayRead(coor, &_coor);
140: c0 = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys);
141: c1 = (xe - 2 - Xs + 1) + (ye - 2 - Ys + 1) * (Xe - Xs) + (ze - 2 - Zs + 1) * (Xe - Xs) * (Ye - Ys);
143: gmin_l[0] = PetscRealPart(_coor[3 * c0 + 0]);
144: gmin_l[1] = PetscRealPart(_coor[3 * c0 + 1]);
145: gmin_l[2] = PetscRealPart(_coor[3 * c0 + 2]);
147: gmax_l[0] = PetscRealPart(_coor[3 * c1 + 0]);
148: gmax_l[1] = PetscRealPart(_coor[3 * c1 + 1]);
149: gmax_l[2] = PetscRealPart(_coor[3 * c1 + 2]);
150: VecRestoreArrayRead(coor, &_coor);
152: if (xs != Xs) xs -= 1;
153: if (ys != Ys) ys -= 1;
154: if (zs != Zs) zs -= 1;
156: mxlocal = xe - xs - 1;
157: mylocal = ye - ys - 1;
158: mzlocal = ze - zs - 1;
160: dx[0] = (gmax_l[0] - gmin_l[0]) / ((PetscReal)mxlocal);
161: dx[1] = (gmax_l[1] - gmin_l[1]) / ((PetscReal)mylocal);
162: dx[2] = (gmax_l[2] - gmin_l[2]) / ((PetscReal)mzlocal);
164: DMGetBoundingBox(dmregular, gmin, gmax);
166: VecGetLocalSize(pos, &n);
167: VecGetBlockSize(pos, &bs);
168: npoints = n / bs;
170: PetscMalloc1(npoints, &cellidx);
171: VecGetArrayRead(pos, &_coor);
172: for (p = 0; p < npoints; p++) {
173: PetscReal coor_p[3];
174: PetscInt mi[3];
176: coor_p[0] = PetscRealPart(_coor[3 * p]);
177: coor_p[1] = PetscRealPart(_coor[3 * p + 1]);
178: coor_p[2] = PetscRealPart(_coor[3 * p + 2]);
180: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
182: if (coor_p[0] < gmin_l[0]) continue;
183: if (coor_p[0] > gmax_l[0]) continue;
184: if (coor_p[1] < gmin_l[1]) continue;
185: if (coor_p[1] > gmax_l[1]) continue;
186: if (coor_p[2] < gmin_l[2]) continue;
187: if (coor_p[2] > gmax_l[2]) continue;
189: for (d = 0; d < 3; d++) mi[d] = (PetscInt)((coor_p[d] - gmin[d]) / dx[d]);
191: if (mi[0] < xs) continue;
192: if (mi[0] > (xe - 1)) continue;
193: if (mi[1] < ys) continue;
194: if (mi[1] > (ye - 1)) continue;
195: if (mi[2] < zs) continue;
196: if (mi[2] > (ze - 1)) continue;
198: if (mi[0] == (xe - 1)) mi[0]--;
199: if (mi[1] == (ye - 1)) mi[1]--;
200: if (mi[2] == (ze - 1)) mi[2]--;
202: cellidx[p] = (mi[0] - xs) + (mi[1] - ys) * mxlocal + (mi[2] - zs) * mxlocal * mylocal;
203: }
204: VecRestoreArrayRead(pos, &_coor);
205: ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell);
206: return 0;
207: }
209: PetscErrorCode DMLocatePoints_DA_Regular(DM dm, Vec pos, DMPointLocationType ltype, PetscSF cellSF)
210: {
211: IS iscell;
212: PetscSFNode *cells;
213: PetscInt p, bs, dim, npoints, nfound;
214: const PetscInt *boxCells;
216: 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: private_DMDALocatePointsIS_2D_Regular(dm, pos, &iscell);
222: break;
223: case 3:
224: private_DMDALocatePointsIS_3D_Regular(dm, pos, &iscell);
225: break;
226: default:
227: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported spatial dimension");
228: }
230: VecGetLocalSize(pos, &npoints);
231: VecGetBlockSize(pos, &bs);
232: npoints = npoints / bs;
234: PetscMalloc1(npoints, &cells);
235: ISGetIndices(iscell, &boxCells);
237: for (p = 0; p < npoints; p++) {
238: cells[p].rank = 0;
239: cells[p].index = boxCells[p];
240: }
241: ISRestoreIndices(iscell, &boxCells);
243: nfound = npoints;
244: PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
245: ISDestroy(&iscell);
246: return 0;
247: }