Actual source code: dageometry.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
5: PETSC_STATIC_INLINE PetscErrorCode FillClosureArray_Private(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, const PetscScalar **array)
6: {
7: PetscScalar *a;
8: PetscInt size = 0, dof, off, d, k, i;
12: for (i = 0; i < nP; ++i) {
13: PetscSectionGetDof(section, points[i], &dof);
14: size += dof;
15: }
16: DMGetWorkArray(dm, size, PETSC_SCALAR, &a);
17: for (i = 0, k = 0; i < nP; ++i) {
18: PetscSectionGetDof(section, points[i], &dof);
19: PetscSectionGetOffset(section, points[i], &off);
21: for (d = 0; d < dof; ++d, ++k) a[k] = vArray[off+d];
22: }
23: *array = a;
24: return(0);
25: }
29: PETSC_STATIC_INLINE PetscErrorCode FillClosureVec_Private(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, const PetscScalar *array, InsertMode mode)
30: {
31: PetscInt dof, off, d, k, i;
35: if ((mode == INSERT_VALUES) || (mode == INSERT_ALL_VALUES)) {
36: for (i = 0, k = 0; i < nP; ++i) {
37: PetscSectionGetDof(section, points[i], &dof);
38: PetscSectionGetOffset(section, points[i], &off);
40: for (d = 0; d < dof; ++d, ++k) vArray[off+d] = array[k];
41: }
42: } else {
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: }
50: return(0);
51: }
55: PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints)
56: {
58: PetscInt *work;
61: if (rn) *rn = n;
62: if (rpoints) {
63: DMGetWorkArray(dm,n,PETSC_INT,&work);
64: PetscMemcpy(work,points,n*sizeof(PetscInt));
65: *rpoints = work;
66: }
67: return(0);
68: }
72: PETSC_STATIC_INLINE PetscErrorCode RestorePointArray_Private(DM dm,PetscInt *rn,const PetscInt **rpoints)
73: {
77: if (rn) *rn = 0;
78: if (rpoints) {
79: DMRestoreWorkArray(dm,*rn, PETSC_INT, (void*) rpoints);
80: }
81: return(0);
82: }
86: PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
87: {
88: DM_DA *da = (DM_DA*) dm->data;
89: PetscInt dim = da->dim;
90: PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
91: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart;
98: if (!section) {DMGetDefaultSection(dm, §ion);}
99: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
100: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
101: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
102: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
103: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
104: DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
105: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
106: xfStart = fStart; xfEnd = xfStart+nXF;
107: yfStart = xfEnd;
108: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
109: if ((p >= cStart) || (p < cEnd)) {
110: /* Cell */
111: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
112: else if (dim == 2) {
113: /* 4 faces, 4 vertices
114: Bottom-left vertex follows same order as cells
115: Bottom y-face same order as cells
116: Left x-face follows same order as cells
117: We number the quad:
119: 8--3--7
120: | |
121: 4 0 2
122: | |
123: 5--1--6
124: */
125: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
126: PetscInt v = cy*nVx + cx + vStart;
127: PetscInt xf = cy*nxF + cx + xfStart;
128: PetscInt yf = c + yfStart;
129: PetscInt points[9];
131: /* Note: initializer list is not computable at compile time, hence we fill the array manually */
132: 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;
134: GetPointArray_Private(dm,9,points,n,closure);
135: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
136: } else if ((p >= vStart) || (p < vEnd)) {
137: /* Vertex */
138: GetPointArray_Private(dm,1,&p,n,closure);
139: } else if ((p >= fStart) || (p < fStart + nXF)) {
140: /* X Face */
141: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
142: else if (dim == 2) {
143: /* 2 vertices: The bottom vertex has the same numbering as the face */
144: PetscInt f = p - xfStart;
145: PetscInt points[3];
147: points[0] = p; points[1] = f; points[2] = f+nVx;
148: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
149: GetPointArray_Private(dm,3,points,n,closure);
150: } else if (dim == 3) {
151: /* 4 vertices */
152: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
153: }
154: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
155: /* Y Face */
156: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
157: else if (dim == 2) {
158: /* 2 vertices: The left vertex has the same numbering as the face */
159: PetscInt f = p - yfStart;
160: PetscInt points[3];
162: points[0] = p; points[1] = f; points[2]= f+1;
163: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
164: GetPointArray_Private(dm, 3, points, n, closure);
165: } else if (dim == 3) {
166: /* 4 vertices */
167: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
168: }
169: } else {
170: /* Z Face */
171: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
172: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
173: else if (dim == 3) {
174: /* 4 vertices */
175: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
176: }
177: }
178: return(0);
179: }
183: PetscErrorCode DMDARestoreClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure)
184: {
191: RestorePointArray_Private(dm,n,closure);
192: return(0);
193: }
197: /* If you did not pass NULL for 'values', you must call DMDARestoreClosureScalar() */
198: PetscErrorCode DMDAGetClosureScalar(DM dm, PetscSection section,PetscInt p,PetscScalar *vArray,const PetscScalar **values)
199: {
200: DM_DA *da = (DM_DA*) dm->data;
201: PetscInt dim = da->dim;
202: PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
203: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;
210: if (!section) {DMGetDefaultSection(dm, §ion);}
211: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
212: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
213: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
214: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
215: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
216: DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
217: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
218: xfStart = fStart; xfEnd = xfStart+nXF;
219: yfStart = xfEnd; yfEnd = yfStart+nYF;
220: zfStart = yfEnd;
221: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
222: if ((p >= cStart) || (p < cEnd)) {
223: /* Cell */
224: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
225: else if (dim == 2) {
226: /* 4 faces, 4 vertices
227: Bottom-left vertex follows same order as cells
228: Bottom y-face same order as cells
229: Left x-face follows same order as cells
230: We number the quad:
232: 8--3--7
233: | |
234: 4 0 2
235: | |
236: 5--1--6
237: */
238: PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1);
239: PetscInt v = cy*nVx + cx + vStart;
240: PetscInt xf = cy*nxF + cx + xfStart;
241: PetscInt yf = c + yfStart;
242: PetscInt points[9];
244: 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;
245: FillClosureArray_Private(dm, section, 9, points, vArray, values);
246: } else {
247: /* 6 faces, 8 vertices
248: Bottom-left-back vertex follows same order as cells
249: Back z-face follows same order as cells
250: Bottom y-face follows same order as cells
251: Left x-face follows same order as cells
253: 14-----13
254: /| /|
255: / | 2 / |
256: / 5| / |
257: 10-----9 4|
258: | 11-|---12
259: |6 / | /
260: | /1 3| /
261: |/ |/
262: 7-----8
263: */
264: PetscInt c = p - cStart;
265: PetscInt points[15];
267: 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;
268: 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;
269: points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;
271: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
272: FillClosureArray_Private(dm, section, 15, points, vArray, values);
273: }
274: } else if ((p >= vStart) || (p < vEnd)) {
275: /* Vertex */
276: FillClosureArray_Private(dm, section, 1, &p, vArray, values);
277: } else if ((p >= fStart) || (p < fStart + nXF)) {
278: /* X 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 bottom vertex has the same numbering as the face */
282: PetscInt f = p - xfStart;
283: PetscInt points[3];
285: points[0] = p; points[1] = f; points[2] = f+nVx;
286: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
287: FillClosureArray_Private(dm, section, 3, points, vArray, values);
288: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
289: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
290: /* Y Face */
291: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
292: else if (dim == 2) {
293: /* 2 vertices: The left vertex has the same numbering as the face */
294: PetscInt f = p - yfStart;
295: PetscInt points[3];
297: points[0] = p; points[1] = f; points[2] = f+1;
298: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken");
299: FillClosureArray_Private(dm, section, 3, points, vArray, values);
300: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
301: } else {
302: /* Z Face */
303: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
304: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
305: else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
306: }
307: return(0);
308: }
312: PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar **values)
313: {
314: PetscScalar *vArray;
321: VecGetArray(v,&vArray);
322: DMDAGetClosureScalar(dm,section,p,vArray,values);
323: VecRestoreArray(v,&vArray);
324: return(0);
325: }
329: PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section,PetscInt p,PetscScalar *vArray,const PetscScalar **values)
330: {
332: PetscInt count;
337: count = 0; /* We are lying about the count */
338: DMRestoreWorkArray(dm, count, PETSC_SCALAR, (void*) values);
339: return(0);
340: }
344: PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar **values)
345: {
352: DMDARestoreClosureScalar(dm,section,p,NULL,values);
353: return(0);
354: }
358: PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode)
359: {
360: DM_DA *da = (DM_DA*) dm->data;
361: PetscInt dim = da->dim;
362: PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF;
363: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart;
370: if (!section) {DMGetDefaultSection(dm, §ion);}
371: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection");
372: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
373: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
374: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
375: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
376: DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);
377: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
378: xfStart = fStart; xfEnd = xfStart+nXF;
379: yfStart = xfEnd; yfEnd = yfStart+nYF;
380: zfStart = yfEnd;
381: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd);
382: if ((p >= cStart) || (p < cEnd)) {
383: /* Cell */
384: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
385: else if (dim == 2) {
386: /* 4 faces, 4 vertices
387: Bottom-left vertex follows same order as cells
388: Bottom y-face same order as cells
389: Left x-face follows same order as cells
390: We number the quad:
392: 8--3--7
393: | |
394: 4 0 2
395: | |
396: 5--1--6
397: */
398: PetscInt c = p - cStart;
399: PetscInt points[9];
401: points[0] = p; points[1] = c+yfStart; points[2] = c+xfStart+1; points[3] = c+yfStart+nyF; points[4] = c+xfStart+0; points[5] = c+vStart+0;
402: points[6] = c+vStart+1; points[7] = c+vStart+nVx+1; points[8] = c+vStart+nVx+0;
403: FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);
404: } else {
405: /* 6 faces, 8 vertices
406: Bottom-left-back vertex follows same order as cells
407: Back z-face follows same order as cells
408: Bottom y-face follows same order as cells
409: Left x-face follows same order as cells
411: 14-----13
412: /| /|
413: / | 2 / |
414: / 5| / |
415: 10-----9 4|
416: | 11-|---12
417: |6 / | /
418: | /1 3| /
419: |/ |/
420: 7-----8
421: */
422: PetscInt c = p - cStart;
423: PetscInt points[15];
425: 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;
426: 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;
427: points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0;
428: FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);
429: }
430: } else if ((p >= vStart) || (p < vEnd)) {
431: /* Vertex */
432: FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);
433: } else if ((p >= fStart) || (p < fStart + nXF)) {
434: /* X Face */
435: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
436: else if (dim == 2) {
437: /* 2 vertices: The bottom vertex has the same numbering as the face */
438: PetscInt f = p - xfStart;
439: PetscInt points[3];
441: points[0] = p; points[1] = f; points[2] = f+nVx;
442: FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
443: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
444: } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) {
445: /* Y Face */
446: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
447: else if (dim == 2) {
448: /* 2 vertices: The left vertex has the same numbering as the face */
449: PetscInt f = p - yfStart;
450: PetscInt points[3];
452: points[0] = p; points[1] = f; points[2] = f+1;
453: FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);
454: } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
455: } else {
456: /* Z Face */
457: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D");
458: else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D");
459: else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented");
460: }
461: return(0);
462: }
466: PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode)
467: {
468: PetscScalar *vArray;
475: VecGetArray(v,&vArray);
476: DMDASetClosureScalar(dm,section,p,vArray,values,mode);
477: VecRestoreArray(v,&vArray);
478: return(0);
479: }
483: /*@
484: DMDAConvertToCell - Convert (i,j,k) to local cell number
486: Not Collective
488: Input Parameter:
489: + da - the distributed array
490: = s - A MatStencil giving (i,j,k)
492: Output Parameter:
493: . cell - the local cell number
495: Level: developer
497: .seealso: DMDAVecGetClosure()
498: @*/
499: PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell)
500: {
501: DM_DA *da = (DM_DA*) dm->data;
502: const PetscInt dim = da->dim;
503: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/;
504: 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;
507: *cell = -1;
508: 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);
509: 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);
510: 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);
511: *cell = (kl*my + jl)*mx + il;
512: return(0);
513: }
517: PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
518: {
519: const PetscScalar x0 = vertices[0];
520: const PetscScalar y0 = vertices[1];
521: const PetscScalar x1 = vertices[2];
522: const PetscScalar y1 = vertices[3];
523: const PetscScalar x2 = vertices[4];
524: const PetscScalar y2 = vertices[5];
525: const PetscScalar x3 = vertices[6];
526: const PetscScalar y3 = vertices[7];
527: const PetscScalar f_01 = x2 - x1 - x3 + x0;
528: const PetscScalar g_01 = y2 - y1 - y3 + y0;
529: const PetscScalar x = refPoint[0];
530: const PetscScalar y = refPoint[1];
531: PetscReal invDet;
532: PetscErrorCode ierr;
535: #if defined(PETSC_USE_DEBUG)
536: PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n",
537: PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));
538: PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));
539: #endif
540: J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5;
541: J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5;
542: *detJ = J[0]*J[3] - J[1]*J[2];
543: invDet = 1.0/(*detJ);
544: invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1];
545: invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0];
546: PetscLogFlops(30);
547: return(0);
548: }
552: PetscErrorCode DMDAComputeCellGeometry(DM dm, PetscInt cell, PetscQuadrature *quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[])
553: {
554: DM cdm;
555: Vec coordinates;
556: const PetscScalar *vertices;
557: PetscInt dim, d, q;
558: PetscErrorCode ierr;
562: DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);
563: DMGetCoordinates(dm, &coordinates);
564: DMGetCoordinateDM(dm, &cdm);
565: DMDAVecGetClosure(cdm, NULL, coordinates, cell, &vertices);
566: for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]);
567: switch (dim) {
568: case 2:
569: for (q = 0; q < quad->numQuadPoints; ++q) {
570: DMDAComputeCellGeometry_2D(dm, vertices, &quad->quadPoints[q*dim], J, invJ, detJ);
571: }
572: break;
573: default:
574: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Dimension %d not supported", dim);
575: }
576: return(0);
577: }