Actual source code: dageometry.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  1: #include <petsc-private/dmdaimpl.h>     /*I  "petscdmda.h"   I*/

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

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

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

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

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

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

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

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

 57:       for (d = 0; d < dof; ++d, ++k) vArray[off+d] += array[k];
 58:     }
 59:   }
 60:   return(0);
 61: }

 65: PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints)
 66: {
 68:   PetscInt       *work;

 71:   if (rn) *rn = n;
 72:   if (rpoints) {
 73:     DMGetWorkArray(dm,n,PETSC_INT,&work);
 74:     PetscMemcpy(work,points,n*sizeof(PetscInt));
 75:     *rpoints = work;
 76:   }
 77:   return(0);
 78: }

 82: PETSC_STATIC_INLINE PetscErrorCode RestorePointArray_Private(DM dm,PetscInt *rn,const PetscInt **rpoints)
 83: {

 87:   if (rn) *rn = 0;
 88:   if (rpoints) {
 89:     DMRestoreWorkArray(dm,*rn, PETSC_INT, (void*) rpoints);
 90:   }
 91:   return(0);
 92: }

 96: PetscErrorCode DMDAGetTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure)
 97: {
 98:   DM_DA         *da = (DM_DA *) dm->data;
 99:   PetscInt       dim = da->dim;
100:   PetscInt       nVx, nVy, nVz, nxF, nXF, nyF, nYF, nzF, nZF;
101:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;

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

126:            8--3--7
127:            |     |
128:            4  0  2
129:            |     |
130:            5--1--6
131:       */
132:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
133:       PetscInt v  = cy*nVx + cx +  vStart;
134:       PetscInt xf = cy*nxF + cx + xfStart;
135:       PetscInt yf = c + yfStart;

138:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
139:       (*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;
140:     } else {
141:       /* 6 faces, 12 edges, 8 vertices
142:          Bottom-left vertex follows same order as cells
143:          Bottom y-face same order as cells
144:          Left x-face follows same order as cells
145:          We number the quad:

147:            8--3--7
148:            |     |
149:            4  0  2
150:            |     |
151:            5--1--6
152:       */
153: #if 0
154:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1), cz = c / ((nVx-1)*(nVy-1));
155:       PetscInt v  = cy*nVx + cx +  vStart;
156:       PetscInt xf = cy*nxF + cx + xfStart;
157:       PetscInt yf = c + yfStart;

160:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
161: #endif
162:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
163:     }
164:   } else if ((p >= vStart) || (p < vEnd)) {
165:     /* Vertex */
167:     if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
168:     (*closure)[0] = p;
169:   } else if ((p >= fStart) || (p < fStart + nXF)) {
170:     /* X Face */
171:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
172:     else if (dim == 2) {
173:       /* 2 vertices: The bottom vertex has the same numbering as the face */
174:       PetscInt f = p - xfStart;

177:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
178:       (*closure)[0] = p; (*closure)[1] = f; (*closure)[2] = f+nVx;
179:     } else if (dim == 3) {
180:       /* 4 vertices */
181:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
182:     }
183:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
184:     /* Y Face */
185:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
186:     else if (dim == 2) {
187:       /* 2 vertices: The left vertex has the same numbering as the face */
188:       PetscInt f = p - yfStart;

191:       if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
192:       (*closure)[0] = p; (*closure)[1] = f; (*closure)[2]= f+1;
193:     } else if (dim == 3) {
194:       /* 4 vertices */
195:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
196:     }
197:   } else {
198:     /* Z Face */
199:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
200:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
201:     else if (dim == 3) {
202:       /* 4 vertices */
203:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
204:     }
205:   }
206:   return(0);
207: }

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

216:   DMRestoreWorkArray(dm, 0, PETSC_INT, closure);
217:   return(0);
218: }

222: PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
223: {
224:   DM_DA          *da = (DM_DA*) dm->data;
225:   PetscInt       dim = da->dim;
226:   PetscInt       nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
227:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;

234:   if (!section) {DMGetDefaultSection(dm, &section);}
235:   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
236:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
237:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
238:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
239:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
240:   DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
241:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
242:   xfStart = fStart; xfEnd = xfStart+nXF;
243:   yfStart = xfEnd;
244:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
245:   if ((p >= cStart) || (p < cEnd)) {
246:     /* Cell */
247:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
248:     else if (dim == 2) {
249:       /* 4 faces, 4 vertices
250:          Bottom-left vertex follows same order as cells
251:          Bottom y-face same order as cells
252:          Left x-face follows same order as cells
253:          We number the quad:

255:            8--3--7
256:            |     |
257:            4  0  2
258:            |     |
259:            5--1--6
260:       */
261:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
262:       PetscInt v  = cy*nVx + cx +  vStart;
263:       PetscInt xf = cy*nxF + cx + xfStart;
264:       PetscInt yf = c + yfStart;
265:       PetscInt points[9];

267:       /* Note: initializer list is not computable at compile time, hence we fill the array manually */
268:       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;

270:       GetPointArray_Private(dm,9,points,n,closure);
271:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
272:   } else if ((p >= vStart) || (p < vEnd)) {
273:     /* Vertex */
274:     GetPointArray_Private(dm,1,&p,n,closure);
275:   } else if ((p >= fStart) || (p < fStart + nXF)) {
276:     /* X Face */
277:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
278:     else if (dim == 2) {
279:       /* 2 vertices: The bottom vertex has the same numbering as the face */
280:       PetscInt f = p - xfStart;
281:       PetscInt points[3];

283:       points[0] = p; points[1] = f; points[2] = f+nVx;
284:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
285:       GetPointArray_Private(dm,3,points,n,closure);
286:     } else if (dim == 3) {
287:       /* 4 vertices */
288:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
289:     }
290:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
291:     /* Y Face */
292:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
293:     else if (dim == 2) {
294:       /* 2 vertices: The left vertex has the same numbering as the face */
295:       PetscInt f = p - yfStart;
296:       PetscInt points[3];

298:       points[0] = p; points[1] = f; points[2]= f+1;
299:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
300:       GetPointArray_Private(dm, 3, points, n, closure);
301:     } else if (dim == 3) {
302:       /* 4 vertices */
303:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
304:     }
305:   } else {
306:     /* Z Face */
307:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
308:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
309:     else if (dim == 3) {
310:       /* 4 vertices */
311:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
312:     }
313:   }
314:   return(0);
315: }

319: /* Zeros n and closure. */
320: PetscErrorCode DMDARestoreClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
321: {

328:   RestorePointArray_Private(dm,n,closure);
329:   return(0);
330: }

334: PetscErrorCode DMDAGetClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values)
335: {
336:   DM_DA          *da = (DM_DA*) dm->data;
337:   PetscInt       dim = da->dim;
338:   PetscInt       nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
339:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;

346:   if (!section) {DMGetDefaultSection(dm, &section);}
347:   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
348:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
349:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
350:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
351:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
352:   DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
353:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
354:   xfStart = fStart; xfEnd = xfStart+nXF;
355:   yfStart = xfEnd;  yfEnd = yfStart+nYF;
356:   zfStart = yfEnd;
357:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
358:   if ((p >= cStart) || (p < cEnd)) {
359:     /* Cell */
360:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
361:     else if (dim == 2) {
362:       /* 4 faces, 4 vertices
363:          Bottom-left vertex follows same order as cells
364:          Bottom y-face same order as cells
365:          Left x-face follows same order as cells
366:          We number the quad:

368:            8--3--7
369:            |     |
370:            4  0  2
371:            |     |
372:            5--1--6
373:       */
374:       PetscInt c  = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
375:       PetscInt v  = cy*nVx + cx +  vStart;
376:       PetscInt xf = cx*nxF + cy + xfStart;
377:       PetscInt yf = c + yfStart;
378:       PetscInt points[9];

380:       points[0] = p;
381:       points[1] = yf;  points[2] = xf+nxF; points[3] = yf+nyF;  points[4] = xf;
382:       points[5] = v+0; points[6] = v+1;    points[7] = v+nVx+1; points[8] = v+nVx+0;
383:       FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);
384:     } else {
385:       /* 6 faces, 8 vertices
386:          Bottom-left-back vertex follows same order as cells
387:          Back z-face follows same order as cells
388:          Bottom y-face follows same order as cells
389:          Left x-face follows same order as cells

391:               14-----13
392:               /|    /|
393:              / | 2 / |
394:             / 5|  /  |
395:           10-----9  4|
396:            |  11-|---12
397:            |6 /  |  /
398:            | /1 3| /
399:            |/    |/
400:            7-----8
401:       */
402:       PetscInt c = p - cStart;
403:       PetscInt points[15];

405:       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;
406:       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;
407:       points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;

409:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
410:       FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);
411:     }
412:   } else if ((p >= vStart) || (p < vEnd)) {
413:     /* Vertex */
414:     FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);
415:   } else if ((p >= fStart) || (p < fStart + nXF)) {
416:     /* X Face */
417:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
418:     else if (dim == 2) {
419:       /* 2 vertices: The bottom vertex has the same numbering as the face */
420:       PetscInt f = p - xfStart;
421:       PetscInt points[3];

423:       points[0] = p; points[1] = f; points[2] = f+nVx;
424:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
425:       FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);
426:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
427:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
428:     /* Y Face */
429:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
430:     else if (dim == 2) {
431:       /* 2 vertices: The left vertex has the same numbering as the face */
432:       PetscInt f = p - yfStart;
433:       PetscInt points[3];

435:       points[0] = p; points[1] = f; points[2] = f+1;
436:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
437:       FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);
438:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
439:   } else {
440:     /* Z Face */
441:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
442:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
443:     else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
444:   }
445:   return(0);
446: }

450: PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values)
451: {
452:   PetscScalar    *vArray;

459:   VecGetArray(v, &vArray);
460:   DMDAGetClosureScalar(dm, section, p, vArray, csize, values);
461:   VecRestoreArray(v, &vArray);
462:   return(0);
463: }

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

474:   DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);
475:   return(0);
476: }

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

488:   DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);
489:   return(0);
490: }

494: PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode)
495: {
496:   DM_DA          *da = (DM_DA*) dm->data;
497:   PetscInt       dim = da->dim;
498:   PetscInt       nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF, nCx, nCy;
499:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;

506:   if (!section) {DMGetDefaultSection(dm, &section);}
507:   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
508:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
509:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
510:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
511:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
512:   DMDAGetNumCells(dm, &nCx, &nCy, NULL, NULL);
513:   DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
514:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
515:   xfStart = fStart; xfEnd = xfStart+nXF;
516:   yfStart = xfEnd;  yfEnd = yfStart+nYF;
517:   zfStart = yfEnd;
518:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
519:   if ((p >= cStart) || (p < cEnd)) {
520:     /* Cell */
521:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
522:     else if (dim == 2) {
523:       /* 4 faces, 4 vertices
524:          Bottom-left vertex follows same order as cells
525:          Bottom y-face same order as cells
526:          Left x-face follows same order as cells
527:          We number the quad:

529:            8--3--7
530:            |     |
531:            4  0  2
532:            |     |
533:            5--1--6
534:       */
535:       PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
536:       PetscInt v  = cy*nVx + cx +  vStart;
537:       PetscInt xf = cx*nxF + cy + xfStart;
538:       PetscInt yf = c + yfStart;
539:       PetscInt points[9];

541:       points[0] = p;
542:       points[1] = yf;  points[2] = xf+nxF; points[3] = yf+nyF;  points[4] = xf;
543:       points[5] = v+0; points[6] = v+1;    points[7] = v+nVx+1; points[8] = v+nVx+0;
544:       FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);
545:     } else {
546:       /* 6 faces, 8 vertices
547:          Bottom-left-back vertex follows same order as cells
548:          Back z-face follows same order as cells
549:          Bottom y-face follows same order as cells
550:          Left x-face follows same order as cells

552:               14-----13
553:               /|    /|
554:              / | 2 / |
555:             / 5|  /  |
556:           10-----9  4|
557:            |  11-|---12
558:            |6 /  |  /
559:            | /1 3| /
560:            |/    |/
561:            7-----8
562:       */
563:       PetscInt c = p - cStart;
564:       PetscInt points[15];

566:       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;
567:       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;
568:       points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;
569:       FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);
570:     }
571:   } else if ((p >= vStart) || (p < vEnd)) {
572:     /* Vertex */
573:     FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);
574:   } else if ((p >= fStart) || (p < fStart + nXF)) {
575:     /* X Face */
576:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
577:     else if (dim == 2) {
578:       /* 2 vertices: The bottom vertex has the same numbering as the face */
579:       PetscInt f = p - xfStart;
580:       PetscInt points[3];

582:       points[0] = p; points[1] = f; points[2] = f+nVx;
583:       FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
584:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
585:   } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
586:     /* Y Face */
587:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
588:     else if (dim == 2) {
589:       /* 2 vertices: The left vertex has the same numbering as the face */
590:       PetscInt f = p - yfStart;
591:       PetscInt points[3];

593:       points[0] = p; points[1] = f; points[2] = f+1;
594:       FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
595:     } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
596:   } else {
597:     /* Z Face */
598:     if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
599:     else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
600:     else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
601:   }
602:   return(0);
603: }

607: PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode)
608: {
609:   PetscScalar    *vArray;

616:   VecGetArray(v,&vArray);
617:   DMDASetClosureScalar(dm,section,p,vArray,values,mode);
618:   VecRestoreArray(v,&vArray);
619:   return(0);
620: }

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

627:   Not Collective

629:   Input Parameter:
630: + da - the distributed array
631: = s - A MatStencil giving (i,j,k)

633:   Output Parameter:
634: . cell - the local cell number

636:   Level: developer

638: .seealso: DMDAVecGetClosure()
639: @*/
640: PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell)
641: {
642:   DM_DA          *da = (DM_DA*) dm->data;
643:   const PetscInt dim = da->dim;
644:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/;
645:   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;

648:   *cell = -1;
649:   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);
650:   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);
651:   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);
652:   *cell = (kl*my + jl)*mx + il;
653:   return(0);
654: }

658: PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
659: {
660:   const PetscScalar x0   = vertices[0];
661:   const PetscScalar y0   = vertices[1];
662:   const PetscScalar x1   = vertices[2];
663:   const PetscScalar y1   = vertices[3];
664:   const PetscScalar x2   = vertices[4];
665:   const PetscScalar y2   = vertices[5];
666:   const PetscScalar x3   = vertices[6];
667:   const PetscScalar y3   = vertices[7];
668:   const PetscScalar f_01 = x2 - x1 - x3 + x0;
669:   const PetscScalar g_01 = y2 - y1 - y3 + y0;
670:   const PetscScalar x    = refPoint[0];
671:   const PetscScalar y    = refPoint[1];
672:   PetscReal         invDet;
673:   PetscErrorCode    ierr;

676: #if 0
677:   PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n",
678:                      PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));
679:   PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));
680: #endif
681:   J[0]    = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5;
682:   J[2]    = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5;
683:   *detJ   = J[0]*J[3] - J[1]*J[2];
684:   invDet  = 1.0/(*detJ);
685:   if (invJ) {
686:     invJ[0] =  invDet*J[3]; invJ[1] = -invDet*J[1];
687:     invJ[2] = -invDet*J[2]; invJ[3] =  invDet*J[0];
688:   }
689:   PetscLogFlops(30);
690:   return(0);
691: }

695: PetscErrorCode DMDAComputeCellGeometry(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[])
696: {
697:   DM               cdm;
698:   Vec              coordinates;
699:   const PetscReal *quadPoints;
700:   PetscScalar     *vertices = NULL;
701:   PetscInt         numQuadPoints, csize, dim, d, q;
702:   PetscErrorCode   ierr;

706:   DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);
707:   DMGetCoordinatesLocal(dm, &coordinates);
708:   DMGetCoordinateDM(dm, &cdm);
709:   DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);
710:   for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]);
711:   switch (dim) {
712:   case 2:
713:     PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, NULL);
714:     for (q = 0; q < numQuadPoints; ++q) {
715:       DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);
716:     }
717:     break;
718:   default:
719:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim);
720:   }
721:   DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);
722:   return(0);
723: }