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