Actual source code: dageometry.c
petsc-3.7.7 2017-09-25
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, §ion);}
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, §ion);}
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, §ion);}
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: }