Actual source code: dageometry.c

petsc-3.6.1 2015-07-22
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:   PetscInt       dim = dm->dim;
 99:   PetscInt       nVx, nVy, nVz, nxF, nXF, nyF, nYF, nzF, nZF;
100:   PetscInt       pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

326:   RestorePointArray_Private(dm,n,closure);
327:   return(0);
328: }

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

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

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

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

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

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

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

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

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

447: PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values)
448: {
449:   PetscScalar    *vArray;

456:   VecGetArray(v, &vArray);
457:   DMDAGetClosureScalar(dm, section, p, vArray, csize, values);
458:   VecRestoreArray(v, &vArray);
459:   return(0);
460: }

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

471:   DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);
472:   return(0);
473: }

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

485:   DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);
486:   return(0);
487: }

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

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

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

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

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

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

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

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

603: PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode)
604: {
605:   PetscScalar    *vArray;

612:   VecGetArray(v,&vArray);
613:   DMDASetClosureScalar(dm,section,p,vArray,values,mode);
614:   VecRestoreArray(v,&vArray);
615:   return(0);
616: }

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

623:   Not Collective

625:   Input Parameter:
626: + da - the distributed array
627: - s - A MatStencil giving (i,j,k)

629:   Output Parameter:
630: . cell - the local cell number

632:   Level: developer

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

644:   *cell = -1;
645:   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);
646:   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);
647:   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);
648:   *cell = (kl*my + jl)*mx + il;
649:   return(0);
650: }

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

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

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

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