Actual source code: dageometry.c
petsc-3.5.4 2015-05-23
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, §ion);}
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, §ion);}
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, §ion);}
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: }