Actual source code: dageometry.c

petsc-3.10.5 2019-03-28
Report Typos and Errors
  1:  #include <petscsf.h>
  2:  #include <petsc/private/dmdaimpl.h>

  4: PETSC_STATIC_INLINE PetscErrorCode FillClosureArray_Static(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, PetscInt *csize, PetscScalar **array)
  5: {
  6:   PetscScalar    *a;
  7:   PetscInt       pStart, pEnd, size = 0, dof, off, d, k, i;

 11:   PetscSectionGetChart(section, &pStart, &pEnd);
 12:   for (i = 0; i < nP; ++i) {
 13:     const PetscInt p = points[i];

 15:     if ((p < pStart) || (p >= pEnd)) continue;
 16:     PetscSectionGetDof(section, p, &dof);
 17:     size += dof;
 18:   }
 19:   if (csize) *csize = size;
 20:   if (array) {
 21:     DMGetWorkArray(dm, size, MPIU_SCALAR, &a);
 22:     for (i = 0, k = 0; i < nP; ++i) {
 23:       const PetscInt p = points[i];

 25:       if ((p < pStart) || (p >= pEnd)) continue;
 26:       PetscSectionGetDof(section, p, &dof);
 27:       PetscSectionGetOffset(section, p, &off);

 29:       for (d = 0; d < dof; ++d, ++k) a[k] = vArray[off+d];
 30:     }
 31:     *array = a;
 32:   }
 33:   return(0);
 34: }

 36: PETSC_STATIC_INLINE PetscErrorCode FillClosureVec_Private(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, const PetscScalar *array, InsertMode mode)
 37: {
 38:   PetscInt       dof, off, d, k, i;

 42:   if ((mode == INSERT_VALUES) || (mode == INSERT_ALL_VALUES)) {
 43:     for (i = 0, k = 0; i < nP; ++i) {
 44:       PetscSectionGetDof(section, points[i], &dof);
 45:       PetscSectionGetOffset(section, points[i], &off);

 47:       for (d = 0; d < dof; ++d, ++k) vArray[off+d] = array[k];
 48:     }
 49:   } else {
 50:     for (i = 0, k = 0; i < nP; ++i) {
 51:       PetscSectionGetDof(section, points[i], &dof);
 52:       PetscSectionGetOffset(section, points[i], &off);

 54:       for (d = 0; d < dof; ++d, ++k) vArray[off+d] += array[k];
 55:     }
 56:   }
 57:   return(0);
 58: }

 60: PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints)
 61: {
 63:   PetscInt       *work;

 66:   if (rn) *rn = n;
 67:   if (rpoints) {
 68:     DMGetWorkArray(dm,n,MPIU_INT,&work);
 69:     PetscMemcpy(work,points,n*sizeof(PetscInt));
 70:     *rpoints = work;
 71:   }
 72:   return(0);
 73: }

 75: PETSC_STATIC_INLINE PetscErrorCode RestorePointArray_Private(DM dm,PetscInt *rn,const PetscInt **rpoints)
 76: {

 80:   if (rn) *rn = 0;
 81:   if (rpoints) {
 82:     /* note the second argument to DMRestoreWorkArray() is always ignored */
 83:     DMRestoreWorkArray(dm,0, MPIU_INT, (void*) rpoints);
 84:   }
 85:   return(0);
 86: }

 88: PetscErrorCode DMDAGetTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure)
 89: {
 90:   PetscInt       dim = dm->dim;
 91:   PetscInt       nVx, nVy, nVz, nxF, nXF, nyF, nYF, nzF, nZF;
 92:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;

 96:   if (!useClosure) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Star operation is not yet supported");
 99:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
100:   DMDAGetHeightStratum(dm,  0,  &cStart, &cEnd);
101:   DMDAGetHeightStratum(dm,  1,  &fStart, &fEnd);
102:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
103:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, NULL);
104:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
105:   xfStart = fStart; xfEnd = xfStart+nXF;
106:   yfStart = xfEnd;
107:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
108:   if ((p >= cStart) || (p < cEnd)) {
109:     /* Cell */
110:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
111:     else if (dim == 2) {
112:       /* 4 faces, 4 vertices
113:          Bottom-left vertex follows same order as cells
114:          Bottom y-face same order as cells
115:          Left x-face follows same order as cells
116:          We number the quad:

118:            8--3--7
119:            |     |
120:            4  0  2
121:            |     |
122:            5--1--6
123:       */
124:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
125:       PetscInt v  = cy*nVx + cx +  vStart;
126:       PetscInt xf = cy*nxF + cx + xfStart;
127:       PetscInt yf = c + yfStart;

129:       *closureSize = 9;
130:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, MPIU_INT, closure);}
131:       (*closure)[0] = p; (*closure)[1] = yf; (*closure)[2] = xf+1; (*closure)[3] = yf+nyF; (*closure)[4] = xf+0; (*closure)[5] = v+0; (*closure)[6]= v+1; (*closure)[7] = v+nVx+1; (*closure)[8] = v+nVx+0;
132:     } else {
133:       /* 6 faces, 12 edges, 8 vertices
134:          Bottom-left vertex follows same order as cells
135:          Bottom y-face same order as cells
136:          Left x-face follows same order as cells
137:          We number the quad:

139:            8--3--7
140:            |     |
141:            4  0  2
142:            |     |
143:            5--1--6
144:       */
145: #if 0
146:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1), cz = c / ((nVx-1)*(nVy-1));
147:       PetscInt v  = cy*nVx + cx +  vStart;
148:       PetscInt xf = cy*nxF + cx + xfStart;
149:       PetscInt yf = c + yfStart;

151:       *closureSize = 26;
152:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, MPIU_INT, closure);}
153: #endif
154:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
155:     }
156:   } else if ((p >= vStart) || (p < vEnd)) {
157:     /* Vertex */
158:     *closureSize = 1;
159:     if (!(*closure)) {DMGetWorkArray(dm, *closureSize, MPIU_INT, closure);}
160:     (*closure)[0] = p;
161:   } else if ((p >= fStart) || (p < fStart + nXF)) {
162:     /* X Face */
163:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
164:     else if (dim == 2) {
165:       /* 2 vertices: The bottom vertex has the same numbering as the face */
166:       PetscInt f = p - xfStart;

168:       *closureSize = 3;
169:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, MPIU_INT, closure);}
170:       (*closure)[0] = p; (*closure)[1] = f; (*closure)[2] = f+nVx;
171:     } else if (dim == 3) {
172:       /* 4 vertices */
173:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
174:     }
175:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
176:     /* Y Face */
177:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
178:     else if (dim == 2) {
179:       /* 2 vertices: The left vertex has the same numbering as the face */
180:       PetscInt f = p - yfStart;

182:       *closureSize = 3;
183:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, MPIU_INT, closure);}
184:       (*closure)[0] = p; (*closure)[1] = f; (*closure)[2]= f+1;
185:     } else if (dim == 3) {
186:       /* 4 vertices */
187:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
188:     }
189:   } else {
190:     /* Z Face */
191:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
192:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
193:     else if (dim == 3) {
194:       /* 4 vertices */
195:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
196:     }
197:   }
198:   return(0);
199: }

201: PetscErrorCode DMDARestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure)
202: {

206:   DMRestoreWorkArray(dm, 0, MPIU_INT, closure);
207:   return(0);
208: }

210: PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
211: {
212:   PetscInt       dim = dm->dim;
213:   PetscInt       nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
214:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;

221:   if (!section) {DMGetSection(dm, &section);}
222:   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
223:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
224:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
225:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
226:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
227:   DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
228:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
229:   xfStart = fStart; xfEnd = xfStart+nXF;
230:   yfStart = xfEnd;
231:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
232:   if ((p >= cStart) || (p < cEnd)) {
233:     /* Cell */
234:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
235:     else if (dim == 2) {
236:       /* 4 faces, 4 vertices
237:          Bottom-left vertex follows same order as cells
238:          Bottom y-face same order as cells
239:          Left x-face follows same order as cells
240:          We number the quad:

242:            8--3--7
243:            |     |
244:            4  0  2
245:            |     |
246:            5--1--6
247:       */
248:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
249:       PetscInt v  = cy*nVx + cx +  vStart;
250:       PetscInt xf = cy*nxF + cx + xfStart;
251:       PetscInt yf = c + yfStart;
252:       PetscInt points[9];

254:       /* Note: initializer list is not computable at compile time, hence we fill the array manually */
255:       points[0] = p; points[1] = yf; points[2] = xf+1; points[3] = yf+nyF; points[4] = xf+0; points[5] = v+0; points[6]= v+1; points[7] = v+nVx+1; points[8] = v+nVx+0;

257:       GetPointArray_Private(dm,9,points,n,closure);
258:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
259:   } else if ((p >= vStart) || (p < vEnd)) {
260:     /* Vertex */
261:     GetPointArray_Private(dm,1,&p,n,closure);
262:   } else if ((p >= fStart) || (p < fStart + nXF)) {
263:     /* X Face */
264:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
265:     else if (dim == 2) {
266:       /* 2 vertices: The bottom vertex has the same numbering as the face */
267:       PetscInt f = p - xfStart;
268:       PetscInt points[3];

270:       points[0] = p; points[1] = f; points[2] = f+nVx;
271:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
272:       GetPointArray_Private(dm,3,points,n,closure);
273:     } else if (dim == 3) {
274:       /* 4 vertices */
275:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
276:     }
277:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
278:     /* Y Face */
279:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
280:     else if (dim == 2) {
281:       /* 2 vertices: The left vertex has the same numbering as the face */
282:       PetscInt f = p - yfStart;
283:       PetscInt points[3];

285:       points[0] = p; points[1] = f; points[2]= f+1;
286:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
287:       GetPointArray_Private(dm, 3, points, n, closure);
288:     } else if (dim == 3) {
289:       /* 4 vertices */
290:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
291:     }
292:   } else {
293:     /* Z Face */
294:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
295:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
296:     else if (dim == 3) {
297:       /* 4 vertices */
298:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
299:     }
300:   }
301:   return(0);
302: }

304: /* Zeros n and closure. */
305: PetscErrorCode DMDARestoreClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
306: {

313:   RestorePointArray_Private(dm,n,closure);
314:   return(0);
315: }

317: PetscErrorCode DMDAGetClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values)
318: {
319:   PetscInt       dim = dm->dim;
320:   PetscInt       nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
321:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;

328:   if (!section) {DMGetSection(dm, &section);}
329:   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
330:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
331:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
332:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
333:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
334:   DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
335:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
336:   xfStart = fStart; xfEnd = xfStart+nXF;
337:   yfStart = xfEnd;  yfEnd = yfStart+nYF;
338:   zfStart = yfEnd;
339:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
340:   if ((p >= cStart) || (p < cEnd)) {
341:     /* Cell */
342:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
343:     else if (dim == 2) {
344:       /* 4 faces, 4 vertices
345:          Bottom-left vertex follows same order as cells
346:          Bottom y-face same order as cells
347:          Left x-face follows same order as cells
348:          We number the quad:

350:            8--3--7
351:            |     |
352:            4  0  2
353:            |     |
354:            5--1--6
355:       */
356:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
357:       PetscInt v  = cy*nVx + cx +  vStart;
358:       PetscInt xf = cx*nxF + cy + xfStart;
359:       PetscInt yf = c + yfStart;
360:       PetscInt points[9];

362:       points[0] = p;
363:       points[1] = yf;  points[2] = xf+nxF; points[3] = yf+nyF;  points[4] = xf;
364:       points[5] = v+0; points[6] = v+1;    points[7] = v+nVx+1; points[8] = v+nVx+0;
365:       FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);
366:     } else {
367:       /* 6 faces, 8 vertices
368:          Bottom-left-back vertex follows same order as cells
369:          Back z-face follows same order as cells
370:          Bottom y-face follows same order as cells
371:          Left x-face follows same order as cells

373:               14-----13
374:               /|    /|
375:              / | 2 / |
376:             / 5|  /  |
377:           10-----9  4|
378:            |  11-|---12
379:            |6 /  |  /
380:            | /1 3| /
381:            |/    |/
382:            7-----8
383:       */
384:       PetscInt c = p - cStart;
385:       PetscInt points[15];

387:       points[0]  = p; points[1] = c+zfStart; points[2] = c+zfStart+nzF; points[3] = c+yfStart; points[4] = c+xfStart+nxF; points[5] = c+yfStart+nyF; points[6] = c+xfStart;
388:       points[7]  = c+vStart+0; points[8] = c+vStart+1; points[9] = c+vStart+nVx+1; points[10] = c+vStart+nVx+0; points[11] = c+vStart+nVx*nVy+0;
389:       points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;

391:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
392:       FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);
393:     }
394:   } else if ((p >= vStart) || (p < vEnd)) {
395:     /* Vertex */
396:     FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);
397:   } else if ((p >= fStart) || (p < fStart + nXF)) {
398:     /* X Face */
399:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
400:     else if (dim == 2) {
401:       /* 2 vertices: The bottom vertex has the same numbering as the face */
402:       PetscInt f = p - xfStart;
403:       PetscInt points[3];

405:       points[0] = p; points[1] = f; points[2] = f+nVx;
406:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
407:       FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);
408:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
409:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
410:     /* Y Face */
411:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
412:     else if (dim == 2) {
413:       /* 2 vertices: The left vertex has the same numbering as the face */
414:       PetscInt f = p - yfStart;
415:       PetscInt points[3];

417:       points[0] = p; points[1] = f; points[2] = f+1;
418:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
419:       FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);
420:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
421:   } else {
422:     /* Z Face */
423:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
424:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
425:     else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
426:   }
427:   return(0);
428: }

430: PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values)
431: {
432:   PetscScalar    *vArray;

439:   VecGetArray(v, &vArray);
440:   DMDAGetClosureScalar(dm, section, p, vArray, csize, values);
441:   VecRestoreArray(v, &vArray);
442:   return(0);
443: }

445: PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values)
446: {

452:   DMRestoreWorkArray(dm, csize ? *csize : 0, MPIU_SCALAR, (void*) values);
453:   return(0);
454: }

456: PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values)
457: {

464:   DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);
465:   return(0);
466: }

468: PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode)
469: {
470:   PetscInt       dim = dm->dim;
471:   PetscInt       nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF, nCx, nCy;
472:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;

479:   if (!section) {DMGetSection(dm, &section);}
480:   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
481:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
482:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
483:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
484:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
485:   DMDAGetNumCells(dm, &nCx, &nCy, NULL, NULL);
486:   DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
487:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
488:   xfStart = fStart; xfEnd = xfStart+nXF;
489:   yfStart = xfEnd;  yfEnd = yfStart+nYF;
490:   zfStart = yfEnd;
491:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
492:   if ((p >= cStart) || (p < cEnd)) {
493:     /* Cell */
494:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
495:     else if (dim == 2) {
496:       /* 4 faces, 4 vertices
497:          Bottom-left vertex follows same order as cells
498:          Bottom y-face same order as cells
499:          Left x-face follows same order as cells
500:          We number the quad:

502:            8--3--7
503:            |     |
504:            4  0  2
505:            |     |
506:            5--1--6
507:       */
508:       PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
509:       PetscInt v  = cy*nVx + cx +  vStart;
510:       PetscInt xf = cx*nxF + cy + xfStart;
511:       PetscInt yf = c + yfStart;
512:       PetscInt points[9];

514:       points[0] = p;
515:       points[1] = yf;  points[2] = xf+nxF; points[3] = yf+nyF;  points[4] = xf;
516:       points[5] = v+0; points[6] = v+1;    points[7] = v+nVx+1; points[8] = v+nVx+0;
517:       FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);
518:     } else {
519:       /* 6 faces, 8 vertices
520:          Bottom-left-back vertex follows same order as cells
521:          Back z-face follows same order as cells
522:          Bottom y-face follows same order as cells
523:          Left x-face follows same order as cells

525:               14-----13
526:               /|    /|
527:              / | 2 / |
528:             / 5|  /  |
529:           10-----9  4|
530:            |  11-|---12
531:            |6 /  |  /
532:            | /1 3| /
533:            |/    |/
534:            7-----8
535:       */
536:       PetscInt c = p - cStart;
537:       PetscInt points[15];

539:       points[0]  = p; points[1] = c+zfStart; points[2] = c+zfStart+nzF; points[3] = c+yfStart; points[4] = c+xfStart+nxF; points[5] = c+yfStart+nyF; points[6] = c+xfStart;
540:       points[7]  = c+vStart+0; points[8] = c+vStart+1; points[9] = c+vStart+nVx+1; points[10] = c+vStart+nVx+0; points[11] = c+vStart+nVx*nVy+0; points[12] = c+vStart+nVx*nVy+1;
541:       points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;
542:       FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);
543:     }
544:   } else if ((p >= vStart) || (p < vEnd)) {
545:     /* Vertex */
546:     FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);
547:   } else if ((p >= fStart) || (p < fStart + nXF)) {
548:     /* X Face */
549:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
550:     else if (dim == 2) {
551:       /* 2 vertices: The bottom vertex has the same numbering as the face */
552:       PetscInt f = p - xfStart;
553:       PetscInt points[3];

555:       points[0] = p; points[1] = f; points[2] = f+nVx;
556:       FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
557:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
558:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
559:     /* Y Face */
560:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
561:     else if (dim == 2) {
562:       /* 2 vertices: The left vertex has the same numbering as the face */
563:       PetscInt f = p - yfStart;
564:       PetscInt points[3];

566:       points[0] = p; points[1] = f; points[2] = f+1;
567:       FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
568:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
569:   } else {
570:     /* Z Face */
571:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
572:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
573:     else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
574:   }
575:   return(0);
576: }

578: PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode)
579: {
580:   PetscScalar    *vArray;

587:   VecGetArray(v,&vArray);
588:   DMDASetClosureScalar(dm,section,p,vArray,values,mode);
589:   VecRestoreArray(v,&vArray);
590:   return(0);
591: }

593: /*@
594:   DMDAConvertToCell - Convert (i,j,k) to local cell number

596:   Not Collective

598:   Input Parameter:
599: + da - the distributed array
600: - s - A MatStencil giving (i,j,k)

602:   Output Parameter:
603: . cell - the local cell number

605:   Level: developer

607: .seealso: DMDAVecGetClosure()
608: @*/
609: PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell)
610: {
611:   DM_DA          *da = (DM_DA*) dm->data;
612:   const PetscInt dim = dm->dim;
613:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/;
614:   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;

617:   *cell = -1;
618:   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);
619:   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);
620:   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);
621:   *cell = (kl*my + jl)*mx + il;
622:   return(0);
623: }

625: PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
626: {
627:   const PetscScalar x0   = vertices[0];
628:   const PetscScalar y0   = vertices[1];
629:   const PetscScalar x1   = vertices[2];
630:   const PetscScalar y1   = vertices[3];
631:   const PetscScalar x2   = vertices[4];
632:   const PetscScalar y2   = vertices[5];
633:   const PetscScalar x3   = vertices[6];
634:   const PetscScalar y3   = vertices[7];
635:   const PetscScalar f_01 = x2 - x1 - x3 + x0;
636:   const PetscScalar g_01 = y2 - y1 - y3 + y0;
637:   const PetscScalar x    = refPoint[0];
638:   const PetscScalar y    = refPoint[1];
639:   PetscReal         invDet;
640:   PetscErrorCode    ierr;

643: #if 0
644:   PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n",
645:                      PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));
646:   PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));
647: #endif
648:   J[0]    = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5;
649:   J[2]    = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5;
650:   *detJ   = J[0]*J[3] - J[1]*J[2];
651:   invDet  = 1.0/(*detJ);
652:   if (invJ) {
653:     invJ[0] =  invDet*J[3]; invJ[1] = -invDet*J[1];
654:     invJ[2] = -invDet*J[2]; invJ[3] =  invDet*J[0];
655:   }
656:   PetscLogFlops(30);
657:   return(0);
658: }

660: PetscErrorCode DMDAComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[])
661: {
662:   DM               cdm;
663:   Vec              coordinates;
664:   const PetscReal *quadPoints;
665:   PetscScalar     *vertices = NULL;
666:   PetscInt         Nq, csize, dim, d, q;
667:   PetscErrorCode   ierr;

671:   DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);
672:   DMGetCoordinatesLocal(dm, &coordinates);
673:   DMGetCoordinateDM(dm, &cdm);
674:   DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);
675:   for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]);
676:   switch (dim) {
677:   case 2:
678:     PetscQuadratureGetData(quad, NULL, NULL, &Nq, &quadPoints, NULL);
679:     for (q = 0; q < Nq; ++q) {
680:       DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);
681:     }
682:     break;
683:   default:
684:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim);
685:   }
686:   DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);
687:   return(0);
688: }

690: PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular,Vec pos,IS *iscell)
691: {
692:   PetscInt          n,bs,p,npoints;
693:   PetscInt          xs,xe,Xs,Xe,mxlocal;
694:   PetscInt          ys,ye,Ys,Ye,mylocal;
695:   PetscInt          d,c0,c1;
696:   PetscReal         gmin_l[2],gmax_l[2],dx[2];
697:   PetscReal         gmin[2],gmax[2];
698:   PetscInt          *cellidx;
699:   Vec               coor;
700:   const PetscScalar *_coor;
701:   PetscErrorCode    ierr;

704:   DMDAGetCorners(dmregular,&xs,&ys,0,&xe,&ye,0);
705:   DMDAGetGhostCorners(dmregular,&Xs,&Ys,0,&Xe,&Ye,0);
706:   xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
707:   ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
708: 
709:   mxlocal = xe - xs - 1;
710:   mylocal = ye - ys - 1;
711: 
712:   DMGetCoordinatesLocal(dmregular,&coor);
713:   VecGetArrayRead(coor,&_coor);
714:   c0 = (xs-Xs) + (ys-Ys)*(Xe-Xs);
715:   c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs);
716: 
717:   gmin_l[0] = PetscRealPart(_coor[2*c0+0]);
718:   gmin_l[1] = PetscRealPart(_coor[2*c0+1]);
719: 
720:   gmax_l[0] = PetscRealPart(_coor[2*c1+0]);
721:   gmax_l[1] = PetscRealPart(_coor[2*c1+1]);
722: 
723:   dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal);
724:   dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal);
725: 
726:   VecRestoreArrayRead(coor,&_coor);
727: 
728:   DMDAGetBoundingBox(dmregular,gmin,gmax);
729: 
730:   VecGetLocalSize(pos,&n);
731:   VecGetBlockSize(pos,&bs);
732:   npoints = n/bs;
733: 
734:   PetscMalloc1(npoints,&cellidx);
735:   VecGetArrayRead(pos,&_coor);
736:   for (p=0; p<npoints; p++) {
737:     PetscReal coor_p[2];
738:     PetscInt  mi[2];
739: 
740:     coor_p[0] = PetscRealPart(_coor[2*p]);
741:     coor_p[1] = PetscRealPart(_coor[2*p+1]);
742: 
743:     cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
744: 
745:     if (coor_p[0] < gmin_l[0]) { continue; }
746:     if (coor_p[0] > gmax_l[0]) { continue; }
747:     if (coor_p[1] < gmin_l[1]) { continue; }
748:     if (coor_p[1] > gmax_l[1]) { continue; }
749: 
750:     for (d=0; d<2; d++) {
751:       mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] );
752:     }
753: 
754:     if (mi[0] < xs)     { continue; }
755:     if (mi[0] > (xe-1)) { continue; }
756:     if (mi[1] < ys)     { continue; }
757:     if (mi[1] > (ye-1)) { continue; }
758: 
759:     if (mi[0] == (xe-1)) { mi[0]--; }
760:     if (mi[1] == (ye-1)) { mi[1]--; }
761: 
762:     cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal;
763:   }
764:   VecRestoreArrayRead(pos,&_coor);
765:   ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
766:   return(0);
767: }

769: PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular,Vec pos,IS *iscell)
770: {
771:   PetscInt          n,bs,p,npoints;
772:   PetscInt          xs,xe,Xs,Xe,mxlocal;
773:   PetscInt          ys,ye,Ys,Ye,mylocal;
774:   PetscInt          zs,ze,Zs,Ze,mzlocal;
775:   PetscInt          d,c0,c1;
776:   PetscReal         gmin_l[3],gmax_l[3],dx[3];
777:   PetscReal         gmin[3],gmax[3];
778:   PetscInt          *cellidx;
779:   Vec               coor;
780:   const PetscScalar *_coor;
781:   PetscErrorCode    ierr;

784:   DMDAGetCorners(dmregular,&xs,&ys,&zs,&xe,&ye,&ze);
785:   DMDAGetGhostCorners(dmregular,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
786:   xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
787:   ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
788:   ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
789: 
790:   mxlocal = xe - xs - 1;
791:   mylocal = ye - ys - 1;
792:   mzlocal = ze - zs - 1;
793: 
794:   DMGetCoordinatesLocal(dmregular,&coor);
795:   VecGetArrayRead(coor,&_coor);
796:   c0 = (xs-Xs)     + (ys-Ys)    *(Xe-Xs) + (zs-Zs)    *(Xe-Xs)*(Ye-Ys);
797:   c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs) + (ze-2-Zs+1)*(Xe-Xs)*(Ye-Ys);
798: 
799:   gmin_l[0] = PetscRealPart(_coor[3*c0+0]);
800:   gmin_l[1] = PetscRealPart(_coor[3*c0+1]);
801:   gmin_l[2] = PetscRealPart(_coor[3*c0+2]);
802: 
803:   gmax_l[0] = PetscRealPart(_coor[3*c1+0]);
804:   gmax_l[1] = PetscRealPart(_coor[3*c1+1]);
805:   gmax_l[2] = PetscRealPart(_coor[3*c1+2]);
806: 
807:   dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal);
808:   dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal);
809:   dx[2] = (gmax_l[2]-gmin_l[2])/((PetscReal)mzlocal);
810: 
811:   VecRestoreArrayRead(coor,&_coor);
812: 
813:   DMDAGetBoundingBox(dmregular,gmin,gmax);
814: 
815:   VecGetLocalSize(pos,&n);
816:   VecGetBlockSize(pos,&bs);
817:   npoints = n/bs;
818: 
819:   PetscMalloc1(npoints,&cellidx);
820:   VecGetArrayRead(pos,&_coor);
821:   for (p=0; p<npoints; p++) {
822:     PetscReal coor_p[3];
823:     PetscInt  mi[3];
824: 
825:     coor_p[0] = PetscRealPart(_coor[3*p]);
826:     coor_p[1] = PetscRealPart(_coor[3*p+1]);
827:     coor_p[2] = PetscRealPart(_coor[3*p+2]);
828: 
829:     cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
830: 
831:     if (coor_p[0] < gmin_l[0]) { continue; }
832:     if (coor_p[0] > gmax_l[0]) { continue; }
833:     if (coor_p[1] < gmin_l[1]) { continue; }
834:     if (coor_p[1] > gmax_l[1]) { continue; }
835:     if (coor_p[2] < gmin_l[2]) { continue; }
836:     if (coor_p[2] > gmax_l[2]) { continue; }
837: 
838:     for (d=0; d<3; d++) {
839:       mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] );
840:     }
841: 
842:     if (mi[0] < xs)     { continue; }
843:     if (mi[0] > (xe-1)) { continue; }
844:     if (mi[1] < ys)     { continue; }
845:     if (mi[1] > (ye-1)) { continue; }
846:     if (mi[2] < zs)     { continue; }
847:     if (mi[2] > (ze-1)) { continue; }
848: 
849:     if (mi[0] == (xe-1)) { mi[0]--; }
850:     if (mi[1] == (ye-1)) { mi[1]--; }
851:     if (mi[2] == (ze-1)) { mi[2]--; }
852: 
853:     cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal + (mi[2]-zs) * mxlocal * mylocal;
854:   }
855:   VecRestoreArrayRead(pos,&_coor);
856:   ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
857:   return(0);
858: }

860: PetscErrorCode DMLocatePoints_DA_Regular(DM dm,Vec pos,DMPointLocationType ltype,PetscSF cellSF)
861: {
862:   IS iscell;
863:   PetscSFNode *cells;
864:   PetscInt p,bs,dim,npoints,nfound;
865:   const PetscInt *boxCells;
867: 
869:   VecGetBlockSize(pos,&dim);
870:   switch (dim) {
871:     case 1:
872:       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support not provided for 1D");
873:       break;
874:     case 2:
875:       private_DMDALocatePointsIS_2D_Regular(dm,pos,&iscell);
876:       break;
877:     case 3:
878:       private_DMDALocatePointsIS_3D_Regular(dm,pos,&iscell);
879:       break;
880:     default:
881:       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupport spatial dimension");
882:       break;
883:   }
884: 
885:   VecGetLocalSize(pos,&npoints);
886:   VecGetBlockSize(pos,&bs);
887:   npoints = npoints / bs;
888: 
889:   PetscMalloc1(npoints, &cells);
890:   ISGetIndices(iscell, &boxCells);
891: 
892:   for (p=0; p<npoints; p++) {
893:     cells[p].rank  = 0;
894:     cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
895: 
896:     cells[p].index = boxCells[p];
897:   }
898:   ISRestoreIndices(iscell, &boxCells);
899: 
900:   nfound = npoints;
901:   PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
902:   ISDestroy(&iscell);
903: 
904:   return(0);
905: }