Actual source code: dageometry.c
petsc-3.13.6 2020-09-29
1: #include <petscsf.h>
2: #include <petsc/private/dmdaimpl.h>
5: /*@
6: DMDAConvertToCell - Convert (i,j,k) to local cell number
8: Not Collective
10: Input Parameter:
11: + da - the distributed array
12: - s - A MatStencil giving (i,j,k)
14: Output Parameter:
15: . cell - the local cell number
17: Level: developer
19: .seealso: DMDAVecGetClosure()
20: @*/
21: PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell)
22: {
23: DM_DA *da = (DM_DA*) dm->data;
24: const PetscInt dim = dm->dim;
25: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/;
26: 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: *cell = -1;
30: if ((s.i < da->Xs/da->w) || (s.i >= da->Xe/da->w)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil i %D should be in [%D, %D)", s.i, da->Xs/da->w, da->Xe/da->w);
31: if ((dim > 1) && ((s.j < da->Ys) || (s.j >= da->Ye))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil j %D should be in [%D, %D)", s.j, da->Ys, da->Ye);
32: if ((dim > 2) && ((s.k < da->Zs) || (s.k >= da->Ze))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil k %D should be in [%D, %D)", s.k, da->Zs, da->Ze);
33: *cell = (kl*my + jl)*mx + il;
34: return(0);
35: }
37: PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular,Vec pos,IS *iscell)
38: {
39: PetscInt n,bs,p,npoints;
40: PetscInt xs,xe,Xs,Xe,mxlocal;
41: PetscInt ys,ye,Ys,Ye,mylocal;
42: PetscInt d,c0,c1;
43: PetscReal gmin_l[2],gmax_l[2],dx[2];
44: PetscReal gmin[2],gmax[2];
45: PetscInt *cellidx;
46: Vec coor;
47: const PetscScalar *_coor;
48: PetscErrorCode ierr;
51: DMDAGetCorners(dmregular,&xs,&ys,0,&xe,&ye,0);
52: DMDAGetGhostCorners(dmregular,&Xs,&Ys,0,&Xe,&Ye,0);
53: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
54: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
56: mxlocal = xe - xs - 1;
57: mylocal = ye - ys - 1;
59: DMGetCoordinatesLocal(dmregular,&coor);
60: VecGetArrayRead(coor,&_coor);
61: c0 = (xs-Xs) + (ys-Ys)*(Xe-Xs);
62: c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs);
64: gmin_l[0] = PetscRealPart(_coor[2*c0+0]);
65: gmin_l[1] = PetscRealPart(_coor[2*c0+1]);
67: gmax_l[0] = PetscRealPart(_coor[2*c1+0]);
68: gmax_l[1] = PetscRealPart(_coor[2*c1+1]);
70: dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal);
71: dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal);
73: VecRestoreArrayRead(coor,&_coor);
75: DMGetBoundingBox(dmregular,gmin,gmax);
77: VecGetLocalSize(pos,&n);
78: VecGetBlockSize(pos,&bs);
79: npoints = n/bs;
81: PetscMalloc1(npoints,&cellidx);
82: VecGetArrayRead(pos,&_coor);
83: for (p=0; p<npoints; p++) {
84: PetscReal coor_p[2];
85: PetscInt mi[2];
87: coor_p[0] = PetscRealPart(_coor[2*p]);
88: coor_p[1] = PetscRealPart(_coor[2*p+1]);
90: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
92: if (coor_p[0] < gmin_l[0]) { continue; }
93: if (coor_p[0] > gmax_l[0]) { continue; }
94: if (coor_p[1] < gmin_l[1]) { continue; }
95: if (coor_p[1] > gmax_l[1]) { continue; }
97: for (d=0; d<2; d++) {
98: mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] );
99: }
101: if (mi[0] < xs) { continue; }
102: if (mi[0] > (xe-1)) { continue; }
103: if (mi[1] < ys) { continue; }
104: if (mi[1] > (ye-1)) { continue; }
106: if (mi[0] == (xe-1)) { mi[0]--; }
107: if (mi[1] == (ye-1)) { mi[1]--; }
109: cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal;
110: }
111: VecRestoreArrayRead(pos,&_coor);
112: ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
113: return(0);
114: }
116: PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular,Vec pos,IS *iscell)
117: {
118: PetscInt n,bs,p,npoints;
119: PetscInt xs,xe,Xs,Xe,mxlocal;
120: PetscInt ys,ye,Ys,Ye,mylocal;
121: PetscInt zs,ze,Zs,Ze,mzlocal;
122: PetscInt d,c0,c1;
123: PetscReal gmin_l[3],gmax_l[3],dx[3];
124: PetscReal gmin[3],gmax[3];
125: PetscInt *cellidx;
126: Vec coor;
127: const PetscScalar *_coor;
128: PetscErrorCode ierr;
131: DMDAGetCorners(dmregular,&xs,&ys,&zs,&xe,&ye,&ze);
132: DMDAGetGhostCorners(dmregular,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
133: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
134: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
135: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
137: mxlocal = xe - xs - 1;
138: mylocal = ye - ys - 1;
139: mzlocal = ze - zs - 1;
141: DMGetCoordinatesLocal(dmregular,&coor);
142: 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]);
154: dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal);
155: dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal);
156: dx[2] = (gmax_l[2]-gmin_l[2])/((PetscReal)mzlocal);
158: VecRestoreArrayRead(coor,&_coor);
160: DMGetBoundingBox(dmregular,gmin,gmax);
162: VecGetLocalSize(pos,&n);
163: VecGetBlockSize(pos,&bs);
164: npoints = n/bs;
166: PetscMalloc1(npoints,&cellidx);
167: VecGetArrayRead(pos,&_coor);
168: for (p=0; p<npoints; p++) {
169: PetscReal coor_p[3];
170: PetscInt mi[3];
172: coor_p[0] = PetscRealPart(_coor[3*p]);
173: coor_p[1] = PetscRealPart(_coor[3*p+1]);
174: coor_p[2] = PetscRealPart(_coor[3*p+2]);
176: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
178: if (coor_p[0] < gmin_l[0]) { continue; }
179: if (coor_p[0] > gmax_l[0]) { continue; }
180: if (coor_p[1] < gmin_l[1]) { continue; }
181: if (coor_p[1] > gmax_l[1]) { continue; }
182: if (coor_p[2] < gmin_l[2]) { continue; }
183: if (coor_p[2] > gmax_l[2]) { continue; }
185: for (d=0; d<3; d++) {
186: mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] );
187: }
189: if (mi[0] < xs) { continue; }
190: if (mi[0] > (xe-1)) { continue; }
191: if (mi[1] < ys) { continue; }
192: if (mi[1] > (ye-1)) { continue; }
193: if (mi[2] < zs) { continue; }
194: if (mi[2] > (ze-1)) { continue; }
196: if (mi[0] == (xe-1)) { mi[0]--; }
197: if (mi[1] == (ye-1)) { mi[1]--; }
198: if (mi[2] == (ze-1)) { mi[2]--; }
200: cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal + (mi[2]-zs) * mxlocal * mylocal;
201: }
202: VecRestoreArrayRead(pos,&_coor);
203: ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
204: return(0);
205: }
207: PetscErrorCode DMLocatePoints_DA_Regular(DM dm,Vec pos,DMPointLocationType ltype,PetscSF cellSF)
208: {
209: IS iscell;
210: PetscSFNode *cells;
211: PetscInt p,bs,dim,npoints,nfound;
212: 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: break;
221: case 2:
222: private_DMDALocatePointsIS_2D_Regular(dm,pos,&iscell);
223: break;
224: case 3:
225: private_DMDALocatePointsIS_3D_Regular(dm,pos,&iscell);
226: break;
227: default:
228: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupport spatial dimension");
229: break;
230: }
232: VecGetLocalSize(pos,&npoints);
233: VecGetBlockSize(pos,&bs);
234: npoints = npoints / bs;
236: PetscMalloc1(npoints, &cells);
237: ISGetIndices(iscell, &boxCells);
239: for (p=0; p<npoints; p++) {
240: cells[p].rank = 0;
241: cells[p].index = boxCells[p];
242: }
243: ISRestoreIndices(iscell, &boxCells);
245: nfound = npoints;
246: PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
247: ISDestroy(&iscell);
248: return(0);
249: }