Actual source code: dageometry.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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: }