Actual source code: dmfieldda.c
1: #include <petsc/private/dmfieldimpl.h>
2: #include <petsc/private/dmimpl.h>
3: #include <petscdmda.h>
5: typedef struct _n_DMField_DA {
6: PetscScalar *cornerVals;
7: PetscScalar *cornerCoeffs;
8: PetscScalar *work;
9: PetscReal coordRange[3][2];
10: } DMField_DA;
12: static PetscErrorCode DMFieldDestroy_DA(DMField field)
13: {
14: DMField_DA *dafield;
16: dafield = (DMField_DA *)field->data;
17: PetscFree3(dafield->cornerVals, dafield->cornerCoeffs, dafield->work);
18: PetscFree(dafield);
19: return 0;
20: }
22: static PetscErrorCode DMFieldView_DA(DMField field, PetscViewer viewer)
23: {
24: DMField_DA *dafield = (DMField_DA *)field->data;
25: PetscBool iascii;
27: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
28: if (iascii) {
29: PetscInt i, c, dim;
30: PetscInt nc;
31: DM dm = field->dm;
33: PetscViewerASCIIPrintf(viewer, "Field corner values:\n");
34: PetscViewerASCIIPushTab(viewer);
35: DMGetDimension(dm, &dim);
36: nc = field->numComponents;
37: for (i = 0, c = 0; i < (1 << dim); i++) {
38: PetscInt j;
40: for (j = 0; j < nc; j++, c++) {
41: PetscScalar val = dafield->cornerVals[nc * i + j];
43: #if !defined(PETSC_USE_COMPLEX)
44: PetscViewerASCIIPrintf(viewer, "%g ", (double)val);
45: #else
46: PetscViewerASCIIPrintf(viewer, "%g+i%g ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val));
47: #endif
48: }
49: PetscViewerASCIIPrintf(viewer, "\n");
50: }
51: PetscViewerASCIIPopTab(viewer);
52: }
53: return 0;
54: }
56: #define MEdot(y, A, x, m, c, cast) \
57: do { \
58: PetscInt _k, _l; \
59: for (_k = 0; _k < (c); _k++) (y)[_k] = 0.; \
60: for (_l = 0; _l < (m); _l++) { \
61: for (_k = 0; _k < (c); _k++) (y)[_k] += cast((A)[(c)*_l + _k]) * (x)[_l]; \
62: } \
63: } while (0)
65: #define MEHess(out, cf, etaB, etaD, dim, nc, cast) \
66: do { \
67: PetscInt _m, _j, _k; \
68: for (_m = 0; _m < (nc) * (dim) * (dim); _m++) (out)[_m] = 0.; \
69: for (_j = 0; _j < (dim); _j++) { \
70: for (_k = _j + 1; _k < (dim); _k++) { \
71: PetscInt _ind = (1 << _j) + (1 << _k); \
72: for (_m = 0; _m < (nc); _m++) { \
73: PetscScalar c = (cf)[_m] * (etaB)[_ind] * (etaD)[_ind]; \
74: (out)[(_m * (dim) + _k) * (dim) + _j] += cast(c); \
75: (out)[(_m * (dim) + _j) * (dim) + _k] += cast(c); \
76: } \
77: } \
78: } \
79: } while (0)
81: static void MultilinearEvaluate(PetscInt dim, PetscReal (*coordRange)[2], PetscInt nc, PetscScalar *cf, PetscScalar *cfWork, PetscInt nPoints, const PetscScalar *points, PetscDataType datatype, void *B, void *D, void *H)
82: {
83: PetscInt i, j, k, l, m;
84: PetscInt whol = 1 << dim;
85: PetscInt half = whol >> 1;
88: if (!B && !D && !H) return;
89: for (i = 0; i < nPoints; i++) {
90: const PetscScalar *point = &points[dim * i];
91: PetscReal deta[3] = {0.};
92: PetscReal etaB[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
93: PetscReal etaD[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
94: PetscReal work[8];
96: for (j = 0; j < dim; j++) {
97: PetscReal e, d;
99: e = (PetscRealPart(point[j]) - coordRange[j][0]) / coordRange[j][1];
100: deta[j] = d = 1. / coordRange[j][1];
101: for (k = 0; k < whol; k++) work[k] = etaB[k];
102: for (k = 0; k < half; k++) {
103: etaB[k] = work[2 * k] * e;
104: etaB[k + half] = work[2 * k + 1];
105: }
106: if (H) {
107: for (k = 0; k < whol; k++) work[k] = etaD[k];
108: for (k = 0; k < half; k++) {
109: etaD[k + half] = work[2 * k];
110: etaD[k] = work[2 * k + 1] * d;
111: }
112: }
113: }
114: if (B) {
115: if (datatype == PETSC_SCALAR) {
116: PetscScalar *out = &((PetscScalar *)B)[nc * i];
118: MEdot(out, cf, etaB, (1 << dim), nc, (PetscScalar));
119: } else {
120: PetscReal *out = &((PetscReal *)B)[nc * i];
122: MEdot(out, cf, etaB, (1 << dim), nc, PetscRealPart);
123: }
124: }
125: if (D) {
126: if (datatype == PETSC_SCALAR) {
127: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
129: for (m = 0; m < nc * dim; m++) out[m] = 0.;
130: } else {
131: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
133: for (m = 0; m < nc * dim; m++) out[m] = 0.;
134: }
135: for (j = 0; j < dim; j++) {
136: PetscReal d = deta[j];
138: for (k = 0; k < whol * nc; k++) cfWork[k] = cf[k];
139: for (k = 0; k < whol; k++) work[k] = etaB[k];
140: for (k = 0; k < half; k++) {
141: PetscReal e;
143: etaB[k] = work[2 * k];
144: etaB[k + half] = e = work[2 * k + 1];
146: for (l = 0; l < nc; l++) {
147: cf[k * nc + l] = cfWork[2 * k * nc + l];
148: cf[(k + half) * nc + l] = cfWork[(2 * k + 1) * nc + l];
149: }
150: if (datatype == PETSC_SCALAR) {
151: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
153: for (l = 0; l < nc; l++) out[l * dim + j] += d * e * cf[k * nc + l];
154: } else {
155: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
157: for (l = 0; l < nc; l++) out[l * dim + j] += d * e * PetscRealPart(cf[k * nc + l]);
158: }
159: }
160: }
161: }
162: if (H) {
163: if (datatype == PETSC_SCALAR) {
164: PetscScalar *out = &((PetscScalar *)H)[nc * dim * dim * i];
166: MEHess(out, cf, etaB, etaD, dim, nc, (PetscScalar));
167: } else {
168: PetscReal *out = &((PetscReal *)H)[nc * dim * dim * i];
170: MEHess(out, cf, etaB, etaD, dim, nc, PetscRealPart);
171: }
172: }
173: }
174: return;
175: }
177: static PetscErrorCode DMFieldEvaluate_DA(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
178: {
179: DM dm;
180: DMField_DA *dafield;
181: PetscInt dim;
182: PetscInt N, n, nc;
183: const PetscScalar *array;
184: PetscReal(*coordRange)[2];
186: dm = field->dm;
187: nc = field->numComponents;
188: dafield = (DMField_DA *)field->data;
189: DMGetDimension(dm, &dim);
190: VecGetLocalSize(points, &N);
192: n = N / dim;
193: coordRange = &(dafield->coordRange[0]);
194: VecGetArrayRead(points, &array);
195: MultilinearEvaluate(dim, coordRange, nc, dafield->cornerCoeffs, dafield->work, n, array, datatype, B, D, H);
196: VecRestoreArrayRead(points, &array);
197: return 0;
198: }
200: static PetscErrorCode DMFieldEvaluateFE_DA(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
201: {
202: PetscInt c, i, j, k, dim, cellsPer[3] = {0}, first[3] = {0}, whol, half;
203: PetscReal stepPer[3] = {0.};
204: PetscReal cellCoordRange[3][2] = {
205: {0., 1.},
206: {0., 1.},
207: {0., 1.}
208: };
209: PetscScalar *cellCoeffs, *work;
210: DM dm;
211: DMDALocalInfo info;
212: PetscInt cStart, cEnd;
213: PetscInt nq, nc;
214: const PetscReal *q;
215: #if defined(PETSC_USE_COMPLEX)
216: PetscScalar *qs;
217: #else
218: const PetscScalar *qs;
219: #endif
220: DMField_DA *dafield;
221: PetscBool isStride;
222: const PetscInt *cells = NULL;
223: PetscInt sfirst = -1, stride = -1, nCells;
225: dafield = (DMField_DA *)field->data;
226: dm = field->dm;
227: nc = field->numComponents;
228: DMDAGetLocalInfo(dm, &info);
229: dim = info.dim;
230: work = dafield->work;
231: stepPer[0] = 1. / info.mx;
232: stepPer[1] = 1. / info.my;
233: stepPer[2] = 1. / info.mz;
234: first[0] = info.gxs;
235: first[1] = info.gys;
236: first[2] = info.gzs;
237: cellsPer[0] = info.gxm;
238: cellsPer[1] = info.gym;
239: cellsPer[2] = info.gzm;
240: /* TODO: probably take components into account */
241: PetscQuadratureGetData(points, NULL, NULL, &nq, &q, NULL);
242: #if defined(PETSC_USE_COMPLEX)
243: DMGetWorkArray(dm, nq * dim, MPIU_SCALAR, &qs);
244: for (i = 0; i < nq * dim; i++) qs[i] = q[i];
245: #else
246: qs = q;
247: #endif
248: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
249: DMGetWorkArray(dm, (1 << dim) * nc, MPIU_SCALAR, &cellCoeffs);
250: whol = (1 << dim);
251: half = whol >> 1;
252: ISGetLocalSize(cellIS, &nCells);
253: PetscObjectTypeCompare((PetscObject)cellIS, ISSTRIDE, &isStride);
254: if (isStride) {
255: ISStrideGetInfo(cellIS, &sfirst, &stride);
256: } else ISGetIndices(cellIS, &cells);
257: for (c = 0; c < nCells; c++) {
258: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
259: PetscInt rem = cell;
260: PetscInt ijk[3] = {0};
261: void *cB, *cD, *cH;
263: if (datatype == PETSC_SCALAR) {
264: cB = B ? &((PetscScalar *)B)[nc * nq * c] : NULL;
265: cD = D ? &((PetscScalar *)D)[nc * nq * dim * c] : NULL;
266: cH = H ? &((PetscScalar *)H)[nc * nq * dim * dim * c] : NULL;
267: } else {
268: cB = B ? &((PetscReal *)B)[nc * nq * c] : NULL;
269: cD = D ? &((PetscReal *)D)[nc * nq * dim * c] : NULL;
270: cH = H ? &((PetscReal *)H)[nc * nq * dim * dim * c] : NULL;
271: }
273: for (i = 0; i < nc * whol; i++) work[i] = dafield->cornerCoeffs[i];
274: for (j = 0; j < dim; j++) {
275: PetscReal e, d;
276: ijk[j] = (rem % cellsPer[j]);
277: rem /= cellsPer[j];
279: e = 2. * (ijk[j] + first[j] + 0.5) * stepPer[j] - 1.;
280: d = stepPer[j];
281: for (i = 0; i < half; i++) {
282: for (k = 0; k < nc; k++) {
283: cellCoeffs[i * nc + k] = work[2 * i * nc + k] * d;
284: cellCoeffs[(i + half) * nc + k] = work[2 * i * nc + k] * e + work[(2 * i + 1) * nc + k];
285: }
286: }
287: for (i = 0; i < whol * nc; i++) work[i] = cellCoeffs[i];
288: }
289: MultilinearEvaluate(dim, cellCoordRange, nc, cellCoeffs, dafield->work, nq, qs, datatype, cB, cD, cH);
290: }
291: if (!isStride) ISRestoreIndices(cellIS, &cells);
292: DMRestoreWorkArray(dm, (1 << dim) * nc, MPIU_SCALAR, &cellCoeffs);
293: #if defined(PETSC_USE_COMPLEX)
294: DMRestoreWorkArray(dm, nq * dim, MPIU_SCALAR, &qs);
295: #endif
296: return 0;
297: }
299: static PetscErrorCode DMFieldEvaluateFV_DA(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
300: {
301: PetscInt c, i, dim, cellsPer[3] = {0}, first[3] = {0};
302: PetscReal stepPer[3] = {0.};
303: DM dm;
304: DMDALocalInfo info;
305: PetscInt cStart, cEnd, numCells;
306: PetscInt nc;
307: PetscScalar *points;
308: DMField_DA *dafield;
309: PetscBool isStride;
310: const PetscInt *cells = NULL;
311: PetscInt sfirst = -1, stride = -1;
313: dafield = (DMField_DA *)field->data;
314: dm = field->dm;
315: nc = field->numComponents;
316: DMDAGetLocalInfo(dm, &info);
317: dim = info.dim;
318: stepPer[0] = 1. / info.mx;
319: stepPer[1] = 1. / info.my;
320: stepPer[2] = 1. / info.mz;
321: first[0] = info.gxs;
322: first[1] = info.gys;
323: first[2] = info.gzs;
324: cellsPer[0] = info.gxm;
325: cellsPer[1] = info.gym;
326: cellsPer[2] = info.gzm;
327: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
328: ISGetLocalSize(cellIS, &numCells);
329: DMGetWorkArray(dm, dim * numCells, MPIU_SCALAR, &points);
330: PetscObjectTypeCompare((PetscObject)cellIS, ISSTRIDE, &isStride);
331: if (isStride) {
332: ISStrideGetInfo(cellIS, &sfirst, &stride);
333: } else ISGetIndices(cellIS, &cells);
334: for (c = 0; c < numCells; c++) {
335: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
336: PetscInt rem = cell;
337: PetscInt ijk[3] = {0};
340: for (i = 0; i < dim; i++) {
341: ijk[i] = (rem % cellsPer[i]);
342: rem /= cellsPer[i];
343: points[dim * c + i] = (ijk[i] + first[i] + 0.5) * stepPer[i];
344: }
345: }
346: if (!isStride) ISRestoreIndices(cellIS, &cells);
347: MultilinearEvaluate(dim, dafield->coordRange, nc, dafield->cornerCoeffs, dafield->work, numCells, points, datatype, B, D, H);
348: DMRestoreWorkArray(dm, dim * numCells, MPIU_SCALAR, &points);
349: return 0;
350: }
352: static PetscErrorCode DMFieldGetDegree_DA(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
353: {
354: DM dm;
355: PetscInt dim, h, imin;
357: dm = field->dm;
358: ISGetMinMax(pointIS, &imin, NULL);
359: DMGetDimension(dm, &dim);
360: for (h = 0; h <= dim; h++) {
361: PetscInt hEnd;
363: DMDAGetHeightStratum(dm, h, NULL, &hEnd);
364: if (imin < hEnd) break;
365: }
366: dim -= h;
367: if (minDegree) *minDegree = 1;
368: if (maxDegree) *maxDegree = dim;
369: return 0;
370: }
372: static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad)
373: {
374: PetscInt h, dim, imax, imin;
375: DM dm;
377: dm = field->dm;
378: ISGetMinMax(cellIS, &imin, &imax);
379: DMGetDimension(dm, &dim);
380: *quad = NULL;
381: for (h = 0; h <= dim; h++) {
382: PetscInt hStart, hEnd;
384: DMDAGetHeightStratum(dm, h, &hStart, &hEnd);
385: if (imin >= hStart && imax < hEnd) break;
386: }
387: dim -= h;
388: if (dim > 0) PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);
390: return 0;
391: }
393: static PetscErrorCode DMFieldInitialize_DA(DMField field)
394: {
395: DM dm;
396: Vec coords = NULL;
397: PetscInt dim, i, j, k;
398: DMField_DA *dafield = (DMField_DA *)field->data;
400: field->ops->destroy = DMFieldDestroy_DA;
401: field->ops->evaluate = DMFieldEvaluate_DA;
402: field->ops->evaluateFE = DMFieldEvaluateFE_DA;
403: field->ops->evaluateFV = DMFieldEvaluateFV_DA;
404: field->ops->getDegree = DMFieldGetDegree_DA;
405: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA;
406: field->ops->view = DMFieldView_DA;
407: dm = field->dm;
408: DMGetDimension(dm, &dim);
409: DMGetCoordinates(dm, &coords);
410: if (coords) {
411: PetscInt n;
412: const PetscScalar *array;
413: PetscReal mins[3][2] = {
414: {PETSC_MAX_REAL, PETSC_MAX_REAL},
415: {PETSC_MAX_REAL, PETSC_MAX_REAL},
416: {PETSC_MAX_REAL, PETSC_MAX_REAL}
417: };
419: VecGetLocalSize(coords, &n);
420: n /= dim;
421: VecGetArrayRead(coords, &array);
422: for (i = 0, k = 0; i < n; i++) {
423: for (j = 0; j < dim; j++, k++) {
424: PetscReal val = PetscRealPart(array[k]);
426: mins[j][0] = PetscMin(mins[j][0], val);
427: mins[j][1] = PetscMin(mins[j][1], -val);
428: }
429: }
430: VecRestoreArrayRead(coords, &array);
431: MPIU_Allreduce((PetscReal *)mins, &(dafield->coordRange[0][0]), 2 * dim, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm));
432: for (j = 0; j < dim; j++) dafield->coordRange[j][1] = -dafield->coordRange[j][1];
433: } else {
434: for (j = 0; j < dim; j++) {
435: dafield->coordRange[j][0] = 0.;
436: dafield->coordRange[j][1] = 1.;
437: }
438: }
439: for (j = 0; j < dim; j++) {
440: PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]);
441: PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]);
443: dafield->coordRange[j][0] = avg;
444: dafield->coordRange[j][1] = dif;
445: }
446: return 0;
447: }
449: PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field)
450: {
451: DMField_DA *dafield;
453: PetscNew(&dafield);
454: field->data = dafield;
455: DMFieldInitialize_DA(field);
456: return 0;
457: }
459: PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues, DMField *field)
460: {
461: DMField b;
462: DMField_DA *dafield;
463: PetscInt dim, nv, i, j, k;
464: PetscInt half;
465: PetscScalar *cv, *cf, *work;
467: DMFieldCreate(dm, nc, DMFIELD_VERTEX, &b);
468: DMFieldSetType(b, DMFIELDDA);
469: dafield = (DMField_DA *)b->data;
470: DMGetDimension(dm, &dim);
471: nv = (1 << dim) * nc;
472: PetscMalloc3(nv, &cv, nv, &cf, nv, &work);
473: for (i = 0; i < nv; i++) cv[i] = cornerValues[i];
474: for (i = 0; i < nv; i++) cf[i] = cv[i];
475: dafield->cornerVals = cv;
476: dafield->cornerCoeffs = cf;
477: dafield->work = work;
478: half = (1 << (dim - 1));
479: for (i = 0; i < dim; i++) {
480: PetscScalar *w;
482: w = work;
483: for (j = 0; j < half; j++) {
484: for (k = 0; k < nc; k++) w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]);
485: }
486: w = &work[j * nc];
487: for (j = 0; j < half; j++) {
488: for (k = 0; k < nc; k++) w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]);
489: }
490: for (j = 0; j < nv; j++) cf[j] = work[j];
491: }
492: *field = b;
493: return 0;
494: }