Actual source code: dageometry.c
petsc-3.8.4 2018-03-24
1: #include <petscsf.h>
2: #include <petsc/private/dmdaimpl.h>
4: PETSC_STATIC_INLINE PetscErrorCode FillClosureArray_Static(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, PetscInt *csize, PetscScalar **array)
5: {
6: PetscScalar *a;
7: PetscInt pStart, pEnd, size = 0, dof, off, d, k, i;
11: PetscSectionGetChart(section, &pStart, &pEnd);
12: for (i = 0; i < nP; ++i) {
13: const PetscInt p = points[i];
15: if ((p < pStart) || (p >= pEnd)) continue;
16: PetscSectionGetDof(section, p, &dof);
17: size += dof;
18: }
19: if (csize) *csize = size;
20: if (array) {
21: DMGetWorkArray(dm, size, PETSC_SCALAR, &a);
22: for (i = 0, k = 0; i < nP; ++i) {
23: const PetscInt p = points[i];
25: if ((p < pStart) || (p >= pEnd)) continue;
26: PetscSectionGetDof(section, p, &dof);
27: PetscSectionGetOffset(section, p, &off);
29: for (d = 0; d < dof; ++d, ++k) a[k] = vArray[off+d];
30: }
31: *array = a;
32: }
33: return(0);
34: }
36: PETSC_STATIC_INLINE PetscErrorCode FillClosureVec_Private(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, const PetscScalar *array, InsertMode mode)
37: {
38: PetscInt dof, off, d, k, i;
42: if ((mode == INSERT_VALUES) || (mode == INSERT_ALL_VALUES)) {
43: for (i = 0, k = 0; i < nP; ++i) {
44: PetscSectionGetDof(section, points[i], &dof);
45: PetscSectionGetOffset(section, points[i], &off);
47: for (d = 0; d < dof; ++d, ++k) vArray[off+d] = array[k];
48: }
49: } else {
50: for (i = 0, k = 0; i < nP; ++i) {
51: PetscSectionGetDof(section, points[i], &dof);
52: PetscSectionGetOffset(section, points[i], &off);
54: for (d = 0; d < dof; ++d, ++k) vArray[off+d] += array[k];
55: }
56: }
57: return(0);
58: }
60: PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints)
61: {
63: PetscInt *work;
66: if (rn) *rn = n;
67: if (rpoints) {
68: DMGetWorkArray(dm,n,PETSC_INT,&work);
69: PetscMemcpy(work,points,n*sizeof(PetscInt));
70: *rpoints = work;
71: }
72: return(0);
73: }
75: PETSC_STATIC_INLINE PetscErrorCode RestorePointArray_Private(DM dm,PetscInt *rn,const PetscInt **rpoints)
76: {
80: if (rn) *rn = 0;
81: if (rpoints) {
82: /* note the second argument to DMRestoreWorkArray() is always ignored */
83: DMRestoreWorkArray(dm,0, PETSC_INT, (void*) rpoints);
84: }
85: return(0);
86: }
88: PetscErrorCode DMDAGetTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure)
89: {
90: PetscInt dim = dm->dim;
91: PetscInt nVx, nVy, nVz, nxF, nXF, nyF, nYF, nzF, nZF;
92: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;
96: if (!useClosure) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Star operation is not yet supported");
99: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
100: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
101: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
102: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
103: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, NULL);
104: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
105: xfStart = fStart; xfEnd = xfStart+nXF;
106: yfStart = xfEnd;
107: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
108: if ((p >= cStart) || (p < cEnd)) {
109: /* Cell */
110: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
111: else if (dim == 2) {
112: /* 4 faces, 4 vertices
113: Bottom-left vertex follows same order as cells
114: Bottom y-face same order as cells
115: Left x-face follows same order as cells
116: We number the quad:
118: 8--3--7
119: | |
120: 4 0 2
121: | |
122: 5--1--6
123: */
124: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
125: PetscInt v = cy*nVx + cx + vStart;
126: PetscInt xf = cy*nxF + cx + xfStart;
127: PetscInt yf = c + yfStart;
129: *closureSize = 9;
130: if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
131: (*closure)[0] = p; (*closure)[1] = yf; (*closure)[2] = xf+1; (*closure)[3] = yf+nyF; (*closure)[4] = xf+0; (*closure)[5] = v+0; (*closure)[6]= v+1; (*closure)[7] = v+nVx+1; (*closure)[8] = v+nVx+0;
132: } else {
133: /* 6 faces, 12 edges, 8 vertices
134: Bottom-left vertex follows same order as cells
135: Bottom y-face same order as cells
136: Left x-face follows same order as cells
137: We number the quad:
139: 8--3--7
140: | |
141: 4 0 2
142: | |
143: 5--1--6
144: */
145: #if 0
146: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1), cz = c / ((nVx-1)*(nVy-1));
147: PetscInt v = cy*nVx + cx + vStart;
148: PetscInt xf = cy*nxF + cx + xfStart;
149: PetscInt yf = c + yfStart;
151: *closureSize = 26;
152: if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
153: #endif
154: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
155: }
156: } else if ((p >= vStart) || (p < vEnd)) {
157: /* Vertex */
158: *closureSize = 1;
159: if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
160: (*closure)[0] = p;
161: } else if ((p >= fStart) || (p < fStart + nXF)) {
162: /* X Face */
163: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
164: else if (dim == 2) {
165: /* 2 vertices: The bottom vertex has the same numbering as the face */
166: PetscInt f = p - xfStart;
168: *closureSize = 3;
169: if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
170: (*closure)[0] = p; (*closure)[1] = f; (*closure)[2] = f+nVx;
171: } else if (dim == 3) {
172: /* 4 vertices */
173: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
174: }
175: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
176: /* Y Face */
177: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
178: else if (dim == 2) {
179: /* 2 vertices: The left vertex has the same numbering as the face */
180: PetscInt f = p - yfStart;
182: *closureSize = 3;
183: if (!(*closure)) {DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);}
184: (*closure)[0] = p; (*closure)[1] = f; (*closure)[2]= f+1;
185: } else if (dim == 3) {
186: /* 4 vertices */
187: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
188: }
189: } else {
190: /* Z Face */
191: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
192: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
193: else if (dim == 3) {
194: /* 4 vertices */
195: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
196: }
197: }
198: return(0);
199: }
201: PetscErrorCode DMDARestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure)
202: {
206: DMRestoreWorkArray(dm, 0, PETSC_INT, closure);
207: return(0);
208: }
210: PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
211: {
212: PetscInt dim = dm->dim;
213: PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
214: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;
221: if (!section) {DMGetDefaultSection(dm, §ion);}
222: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
223: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
224: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
225: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
226: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
227: DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
228: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
229: xfStart = fStart; xfEnd = xfStart+nXF;
230: yfStart = xfEnd;
231: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
232: if ((p >= cStart) || (p < cEnd)) {
233: /* Cell */
234: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
235: else if (dim == 2) {
236: /* 4 faces, 4 vertices
237: Bottom-left vertex follows same order as cells
238: Bottom y-face same order as cells
239: Left x-face follows same order as cells
240: We number the quad:
242: 8--3--7
243: | |
244: 4 0 2
245: | |
246: 5--1--6
247: */
248: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
249: PetscInt v = cy*nVx + cx + vStart;
250: PetscInt xf = cy*nxF + cx + xfStart;
251: PetscInt yf = c + yfStart;
252: PetscInt points[9];
254: /* Note: initializer list is not computable at compile time, hence we fill the array manually */
255: points[0] = p; points[1] = yf; points[2] = xf+1; points[3] = yf+nyF; points[4] = xf+0; points[5] = v+0; points[6]= v+1; points[7] = v+nVx+1; points[8] = v+nVx+0;
257: GetPointArray_Private(dm,9,points,n,closure);
258: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
259: } else if ((p >= vStart) || (p < vEnd)) {
260: /* Vertex */
261: GetPointArray_Private(dm,1,&p,n,closure);
262: } else if ((p >= fStart) || (p < fStart + nXF)) {
263: /* X Face */
264: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
265: else if (dim == 2) {
266: /* 2 vertices: The bottom vertex has the same numbering as the face */
267: PetscInt f = p - xfStart;
268: PetscInt points[3];
270: points[0] = p; points[1] = f; points[2] = f+nVx;
271: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
272: GetPointArray_Private(dm,3,points,n,closure);
273: } else if (dim == 3) {
274: /* 4 vertices */
275: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
276: }
277: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
278: /* Y Face */
279: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
280: else if (dim == 2) {
281: /* 2 vertices: The left vertex has the same numbering as the face */
282: PetscInt f = p - yfStart;
283: PetscInt points[3];
285: points[0] = p; points[1] = f; points[2]= f+1;
286: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
287: GetPointArray_Private(dm, 3, points, n, closure);
288: } else if (dim == 3) {
289: /* 4 vertices */
290: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
291: }
292: } else {
293: /* Z Face */
294: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
295: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
296: else if (dim == 3) {
297: /* 4 vertices */
298: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
299: }
300: }
301: return(0);
302: }
304: /* Zeros n and closure. */
305: PetscErrorCode DMDARestoreClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
306: {
313: RestorePointArray_Private(dm,n,closure);
314: return(0);
315: }
317: PetscErrorCode DMDAGetClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values)
318: {
319: PetscInt dim = dm->dim;
320: PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
321: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;
328: if (!section) {DMGetDefaultSection(dm, §ion);}
329: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
330: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
331: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
332: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
333: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
334: DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
335: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
336: xfStart = fStart; xfEnd = xfStart+nXF;
337: yfStart = xfEnd; yfEnd = yfStart+nYF;
338: zfStart = yfEnd;
339: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
340: if ((p >= cStart) || (p < cEnd)) {
341: /* Cell */
342: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
343: else if (dim == 2) {
344: /* 4 faces, 4 vertices
345: Bottom-left vertex follows same order as cells
346: Bottom y-face same order as cells
347: Left x-face follows same order as cells
348: We number the quad:
350: 8--3--7
351: | |
352: 4 0 2
353: | |
354: 5--1--6
355: */
356: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
357: PetscInt v = cy*nVx + cx + vStart;
358: PetscInt xf = cx*nxF + cy + xfStart;
359: PetscInt yf = c + yfStart;
360: PetscInt points[9];
362: points[0] = p;
363: points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf;
364: points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0;
365: FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);
366: } else {
367: /* 6 faces, 8 vertices
368: Bottom-left-back vertex follows same order as cells
369: Back z-face follows same order as cells
370: Bottom y-face follows same order as cells
371: Left x-face follows same order as cells
373: 14-----13
374: /| /|
375: / | 2 / |
376: / 5| / |
377: 10-----9 4|
378: | 11-|---12
379: |6 / | /
380: | /1 3| /
381: |/ |/
382: 7-----8
383: */
384: PetscInt c = p - cStart;
385: PetscInt points[15];
387: points[0] = p; points[1] = c+zfStart; points[2] = c+zfStart+nzF; points[3] = c+yfStart; points[4] = c+xfStart+nxF; points[5] = c+yfStart+nyF; points[6] = c+xfStart;
388: points[7] = c+vStart+0; points[8] = c+vStart+1; points[9] = c+vStart+nVx+1; points[10] = c+vStart+nVx+0; points[11] = c+vStart+nVx*nVy+0;
389: points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;
391: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
392: FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);
393: }
394: } else if ((p >= vStart) || (p < vEnd)) {
395: /* Vertex */
396: FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);
397: } else if ((p >= fStart) || (p < fStart + nXF)) {
398: /* X Face */
399: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
400: else if (dim == 2) {
401: /* 2 vertices: The bottom vertex has the same numbering as the face */
402: PetscInt f = p - xfStart;
403: PetscInt points[3];
405: points[0] = p; points[1] = f; points[2] = f+nVx;
406: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
407: FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);
408: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
409: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
410: /* Y Face */
411: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
412: else if (dim == 2) {
413: /* 2 vertices: The left vertex has the same numbering as the face */
414: PetscInt f = p - yfStart;
415: PetscInt points[3];
417: points[0] = p; points[1] = f; points[2] = f+1;
418: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
419: FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);
420: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
421: } else {
422: /* Z Face */
423: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
424: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
425: else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
426: }
427: return(0);
428: }
430: PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values)
431: {
432: PetscScalar *vArray;
439: VecGetArray(v, &vArray);
440: DMDAGetClosureScalar(dm, section, p, vArray, csize, values);
441: VecRestoreArray(v, &vArray);
442: return(0);
443: }
445: PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values)
446: {
452: DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);
453: return(0);
454: }
456: PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values)
457: {
464: DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);
465: return(0);
466: }
468: PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode)
469: {
470: PetscInt dim = dm->dim;
471: PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF, nCx, nCy;
472: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;
479: if (!section) {DMGetDefaultSection(dm, §ion);}
480: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
481: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
482: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
483: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
484: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
485: DMDAGetNumCells(dm, &nCx, &nCy, NULL, NULL);
486: DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
487: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
488: xfStart = fStart; xfEnd = xfStart+nXF;
489: yfStart = xfEnd; yfEnd = yfStart+nYF;
490: zfStart = yfEnd;
491: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
492: if ((p >= cStart) || (p < cEnd)) {
493: /* Cell */
494: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
495: else if (dim == 2) {
496: /* 4 faces, 4 vertices
497: Bottom-left vertex follows same order as cells
498: Bottom y-face same order as cells
499: Left x-face follows same order as cells
500: We number the quad:
502: 8--3--7
503: | |
504: 4 0 2
505: | |
506: 5--1--6
507: */
508: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
509: PetscInt v = cy*nVx + cx + vStart;
510: PetscInt xf = cx*nxF + cy + xfStart;
511: PetscInt yf = c + yfStart;
512: PetscInt points[9];
514: points[0] = p;
515: points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf;
516: points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0;
517: FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);
518: } else {
519: /* 6 faces, 8 vertices
520: Bottom-left-back vertex follows same order as cells
521: Back z-face follows same order as cells
522: Bottom y-face follows same order as cells
523: Left x-face follows same order as cells
525: 14-----13
526: /| /|
527: / | 2 / |
528: / 5| / |
529: 10-----9 4|
530: | 11-|---12
531: |6 / | /
532: | /1 3| /
533: |/ |/
534: 7-----8
535: */
536: PetscInt c = p - cStart;
537: PetscInt points[15];
539: points[0] = p; points[1] = c+zfStart; points[2] = c+zfStart+nzF; points[3] = c+yfStart; points[4] = c+xfStart+nxF; points[5] = c+yfStart+nyF; points[6] = c+xfStart;
540: points[7] = c+vStart+0; points[8] = c+vStart+1; points[9] = c+vStart+nVx+1; points[10] = c+vStart+nVx+0; points[11] = c+vStart+nVx*nVy+0; points[12] = c+vStart+nVx*nVy+1;
541: points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;
542: FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);
543: }
544: } else if ((p >= vStart) || (p < vEnd)) {
545: /* Vertex */
546: FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);
547: } else if ((p >= fStart) || (p < fStart + nXF)) {
548: /* X Face */
549: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
550: else if (dim == 2) {
551: /* 2 vertices: The bottom vertex has the same numbering as the face */
552: PetscInt f = p - xfStart;
553: PetscInt points[3];
555: points[0] = p; points[1] = f; points[2] = f+nVx;
556: FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
557: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
558: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
559: /* Y Face */
560: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
561: else if (dim == 2) {
562: /* 2 vertices: The left vertex has the same numbering as the face */
563: PetscInt f = p - yfStart;
564: PetscInt points[3];
566: points[0] = p; points[1] = f; points[2] = f+1;
567: FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
568: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
569: } else {
570: /* Z Face */
571: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
572: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
573: else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
574: }
575: return(0);
576: }
578: PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode)
579: {
580: PetscScalar *vArray;
587: VecGetArray(v,&vArray);
588: DMDASetClosureScalar(dm,section,p,vArray,values,mode);
589: VecRestoreArray(v,&vArray);
590: return(0);
591: }
593: /*@
594: DMDAConvertToCell - Convert (i,j,k) to local cell number
596: Not Collective
598: Input Parameter:
599: + da - the distributed array
600: - s - A MatStencil giving (i,j,k)
602: Output Parameter:
603: . cell - the local cell number
605: Level: developer
607: .seealso: DMDAVecGetClosure()
608: @*/
609: PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell)
610: {
611: DM_DA *da = (DM_DA*) dm->data;
612: const PetscInt dim = dm->dim;
613: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/;
614: const PetscInt il = s.i - da->Xs/da->w, jl = dim > 1 ? s.j - da->Ys : 0, kl = dim > 2 ? s.k - da->Zs : 0;
617: *cell = -1;
618: if ((s.i < da->Xs/da->w) || (s.i >= da->Xe/da->w)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil i %d should be in [%d, %d)", s.i, da->Xs/da->w, da->Xe/da->w);
619: if ((dim > 1) && ((s.j < da->Ys) || (s.j >= da->Ye))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil j %d should be in [%d, %d)", s.j, da->Ys, da->Ye);
620: if ((dim > 2) && ((s.k < da->Zs) || (s.k >= da->Ze))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil k %d should be in [%d, %d)", s.k, da->Zs, da->Ze);
621: *cell = (kl*my + jl)*mx + il;
622: return(0);
623: }
625: PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
626: {
627: const PetscScalar x0 = vertices[0];
628: const PetscScalar y0 = vertices[1];
629: const PetscScalar x1 = vertices[2];
630: const PetscScalar y1 = vertices[3];
631: const PetscScalar x2 = vertices[4];
632: const PetscScalar y2 = vertices[5];
633: const PetscScalar x3 = vertices[6];
634: const PetscScalar y3 = vertices[7];
635: const PetscScalar f_01 = x2 - x1 - x3 + x0;
636: const PetscScalar g_01 = y2 - y1 - y3 + y0;
637: const PetscScalar x = refPoint[0];
638: const PetscScalar y = refPoint[1];
639: PetscReal invDet;
640: PetscErrorCode ierr;
643: #if 0
644: PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n",
645: PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));
646: PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));
647: #endif
648: J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5;
649: J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5;
650: *detJ = J[0]*J[3] - J[1]*J[2];
651: invDet = 1.0/(*detJ);
652: if (invJ) {
653: invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1];
654: invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0];
655: }
656: PetscLogFlops(30);
657: return(0);
658: }
660: PetscErrorCode DMDAComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[])
661: {
662: DM cdm;
663: Vec coordinates;
664: const PetscReal *quadPoints;
665: PetscScalar *vertices = NULL;
666: PetscInt Nq, csize, dim, d, q;
667: PetscErrorCode ierr;
671: DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);
672: DMGetCoordinatesLocal(dm, &coordinates);
673: DMGetCoordinateDM(dm, &cdm);
674: DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);
675: for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]);
676: switch (dim) {
677: case 2:
678: PetscQuadratureGetData(quad, NULL, NULL, &Nq, &quadPoints, NULL);
679: for (q = 0; q < Nq; ++q) {
680: DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);
681: }
682: break;
683: default:
684: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim);
685: }
686: DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);
687: return(0);
688: }
690: PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular,Vec pos,IS *iscell)
691: {
692: PetscInt n,bs,p,npoints;
693: PetscInt xs,xe,Xs,Xe,mxlocal;
694: PetscInt ys,ye,Ys,Ye,mylocal;
695: PetscInt d,c0,c1;
696: PetscReal gmin_l[2],gmax_l[2],dx[2];
697: PetscReal gmin[2],gmax[2];
698: PetscInt *cellidx;
699: Vec coor;
700: const PetscScalar *_coor;
701: PetscErrorCode ierr;
704: DMDAGetCorners(dmregular,&xs,&ys,0,&xe,&ye,0);
705: DMDAGetGhostCorners(dmregular,&Xs,&Ys,0,&Xe,&Ye,0);
706: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
707: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
708:
709: mxlocal = xe - xs - 1;
710: mylocal = ye - ys - 1;
711:
712: DMGetCoordinatesLocal(dmregular,&coor);
713: VecGetArrayRead(coor,&_coor);
714: c0 = (xs-Xs) + (ys-Ys)*(Xe-Xs);
715: c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs);
716:
717: gmin_l[0] = PetscRealPart(_coor[2*c0+0]);
718: gmin_l[1] = PetscRealPart(_coor[2*c0+1]);
719:
720: gmax_l[0] = PetscRealPart(_coor[2*c1+0]);
721: gmax_l[1] = PetscRealPart(_coor[2*c1+1]);
722:
723: dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal);
724: dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal);
725:
726: VecRestoreArrayRead(coor,&_coor);
727:
728: DMDAGetBoundingBox(dmregular,gmin,gmax);
729:
730: VecGetLocalSize(pos,&n);
731: VecGetBlockSize(pos,&bs);
732: npoints = n/bs;
733:
734: PetscMalloc1(npoints,&cellidx);
735: VecGetArrayRead(pos,&_coor);
736: for (p=0; p<npoints; p++) {
737: PetscReal coor_p[2];
738: PetscInt mi[2];
739:
740: coor_p[0] = PetscRealPart(_coor[2*p]);
741: coor_p[1] = PetscRealPart(_coor[2*p+1]);
742:
743: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
744:
745: if (coor_p[0] < gmin_l[0]) { continue; }
746: if (coor_p[0] > gmax_l[0]) { continue; }
747: if (coor_p[1] < gmin_l[1]) { continue; }
748: if (coor_p[1] > gmax_l[1]) { continue; }
749:
750: for (d=0; d<2; d++) {
751: mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] );
752: }
753:
754: if (mi[0] < xs) { continue; }
755: if (mi[0] > (xe-1)) { continue; }
756: if (mi[1] < ys) { continue; }
757: if (mi[1] > (ye-1)) { continue; }
758:
759: if (mi[0] == (xe-1)) { mi[0]--; }
760: if (mi[1] == (ye-1)) { mi[1]--; }
761:
762: cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal;
763: }
764: VecRestoreArrayRead(pos,&_coor);
765: ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
766: return(0);
767: }
769: PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular,Vec pos,IS *iscell)
770: {
771: PetscInt n,bs,p,npoints;
772: PetscInt xs,xe,Xs,Xe,mxlocal;
773: PetscInt ys,ye,Ys,Ye,mylocal;
774: PetscInt zs,ze,Zs,Ze,mzlocal;
775: PetscInt d,c0,c1;
776: PetscReal gmin_l[3],gmax_l[3],dx[3];
777: PetscReal gmin[3],gmax[3];
778: PetscInt *cellidx;
779: Vec coor;
780: const PetscScalar *_coor;
781: PetscErrorCode ierr;
784: DMDAGetCorners(dmregular,&xs,&ys,&zs,&xe,&ye,&ze);
785: DMDAGetGhostCorners(dmregular,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
786: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
787: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
788: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
789:
790: mxlocal = xe - xs - 1;
791: mylocal = ye - ys - 1;
792: mzlocal = ze - zs - 1;
793:
794: DMGetCoordinatesLocal(dmregular,&coor);
795: VecGetArrayRead(coor,&_coor);
796: c0 = (xs-Xs) + (ys-Ys) *(Xe-Xs) + (zs-Zs) *(Xe-Xs)*(Ye-Ys);
797: c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs) + (ze-2-Zs+1)*(Xe-Xs)*(Ye-Ys);
798:
799: gmin_l[0] = PetscRealPart(_coor[3*c0+0]);
800: gmin_l[1] = PetscRealPart(_coor[3*c0+1]);
801: gmin_l[2] = PetscRealPart(_coor[3*c0+2]);
802:
803: gmax_l[0] = PetscRealPart(_coor[3*c1+0]);
804: gmax_l[1] = PetscRealPart(_coor[3*c1+1]);
805: gmax_l[2] = PetscRealPart(_coor[3*c1+2]);
806:
807: dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal);
808: dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal);
809: dx[2] = (gmax_l[2]-gmin_l[2])/((PetscReal)mzlocal);
810:
811: VecRestoreArrayRead(coor,&_coor);
812:
813: DMDAGetBoundingBox(dmregular,gmin,gmax);
814:
815: VecGetLocalSize(pos,&n);
816: VecGetBlockSize(pos,&bs);
817: npoints = n/bs;
818:
819: PetscMalloc1(npoints,&cellidx);
820: VecGetArrayRead(pos,&_coor);
821: for (p=0; p<npoints; p++) {
822: PetscReal coor_p[3];
823: PetscInt mi[3];
824:
825: coor_p[0] = PetscRealPart(_coor[3*p]);
826: coor_p[1] = PetscRealPart(_coor[3*p+1]);
827: coor_p[2] = PetscRealPart(_coor[3*p+2]);
828:
829: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
830:
831: if (coor_p[0] < gmin_l[0]) { continue; }
832: if (coor_p[0] > gmax_l[0]) { continue; }
833: if (coor_p[1] < gmin_l[1]) { continue; }
834: if (coor_p[1] > gmax_l[1]) { continue; }
835: if (coor_p[2] < gmin_l[2]) { continue; }
836: if (coor_p[2] > gmax_l[2]) { continue; }
837:
838: for (d=0; d<3; d++) {
839: mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] );
840: }
841:
842: if (mi[0] < xs) { continue; }
843: if (mi[0] > (xe-1)) { continue; }
844: if (mi[1] < ys) { continue; }
845: if (mi[1] > (ye-1)) { continue; }
846: if (mi[2] < zs) { continue; }
847: if (mi[2] > (ze-1)) { continue; }
848:
849: if (mi[0] == (xe-1)) { mi[0]--; }
850: if (mi[1] == (ye-1)) { mi[1]--; }
851: if (mi[2] == (ze-1)) { mi[2]--; }
852:
853: cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal + (mi[2]-zs) * mxlocal * mylocal;
854: }
855: VecRestoreArrayRead(pos,&_coor);
856: ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
857: return(0);
858: }
860: PetscErrorCode DMLocatePoints_DA_Regular(DM dm,Vec pos,DMPointLocationType ltype,PetscSF cellSF)
861: {
862: IS iscell;
863: PetscSFNode *cells;
864: PetscInt p,bs,dim,npoints,nfound;
865: const PetscInt *boxCells;
867:
869: VecGetBlockSize(pos,&dim);
870: switch (dim) {
871: case 1:
872: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support not provided for 1D");
873: break;
874: case 2:
875: private_DMDALocatePointsIS_2D_Regular(dm,pos,&iscell);
876: break;
877: case 3:
878: private_DMDALocatePointsIS_3D_Regular(dm,pos,&iscell);
879: break;
880: default:
881: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupport spatial dimension");
882: break;
883: }
884:
885: VecGetLocalSize(pos,&npoints);
886: VecGetBlockSize(pos,&bs);
887: npoints = npoints / bs;
888:
889: PetscMalloc1(npoints, &cells);
890: ISGetIndices(iscell, &boxCells);
891:
892: for (p=0; p<npoints; p++) {
893: cells[p].rank = 0;
894: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
895:
896: cells[p].index = boxCells[p];
897: }
898: ISRestoreIndices(iscell, &boxCells);
899:
900: nfound = npoints;
901: PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
902: ISDestroy(&iscell);
903:
904: return(0);
905: }