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