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: }