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: PetscFunctionBegin;
17: dafield = (DMField_DA *)field->data;
18: PetscCall(PetscFree3(dafield->cornerVals, dafield->cornerCoeffs, dafield->work));
19: PetscCall(PetscFree(dafield));
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: static PetscErrorCode DMFieldView_DA(DMField field, PetscViewer viewer)
24: {
25: DMField_DA *dafield = (DMField_DA *)field->data;
26: PetscBool isascii;
28: PetscFunctionBegin;
29: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
30: if (isascii) {
31: PetscInt i, dim;
32: PetscInt nc;
33: DM dm = field->dm;
35: PetscCall(PetscViewerASCIIPrintf(viewer, "Field corner values:\n"));
36: PetscCall(PetscViewerASCIIPushTab(viewer));
37: PetscCall(DMGetDimension(dm, &dim));
38: nc = field->numComponents;
39: for (i = 0; i < (1 << dim); i++) {
40: PetscInt j;
42: for (j = 0; j < nc; j++) {
43: PetscScalar val = dafield->cornerVals[nc * i + j];
45: #if !PetscDefined(USE_COMPLEX)
46: PetscCall(PetscViewerASCIIPrintf(viewer, "%g ", (double)val));
47: #else
48: PetscCall(PetscViewerASCIIPrintf(viewer, "%g+i%g ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
49: #endif
50: }
51: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
52: }
53: PetscCall(PetscViewerASCIIPopTab(viewer));
54: }
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: #define MEdot(y, A, x, m, c, cast) \
59: do { \
60: PetscInt _k, _l; \
61: for (_k = 0; _k < (c); _k++) (y)[_k] = 0.; \
62: for (_l = 0; _l < (m); _l++) { \
63: for (_k = 0; _k < (c); _k++) (y)[_k] += cast((A)[(c) * _l + _k]) * (x)[_l]; \
64: } \
65: } while (0)
67: #define MEHess(out, cf, etaB, etaD, dim, nc, cast) \
68: do { \
69: PetscInt _m, _j, _k; \
70: for (_m = 0; _m < (nc) * (dim) * (dim); _m++) (out)[_m] = 0.; \
71: for (_j = 0; _j < (dim); _j++) { \
72: for (_k = _j + 1; _k < (dim); _k++) { \
73: PetscInt _ind = (1 << _j) + (1 << _k); \
74: for (_m = 0; _m < (nc); _m++) { \
75: PetscScalar c = (cf)[_m] * (etaB)[_ind] * (etaD)[_ind]; \
76: (out)[(_m * (dim) + _k) * (dim) + _j] += cast(c); \
77: (out)[(_m * (dim) + _j) * (dim) + _k] += cast(c); \
78: } \
79: } \
80: } \
81: } while (0)
83: 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)
84: {
85: PetscInt i, j, k, l, m;
86: PetscInt whol = 1 << dim;
87: PetscInt half = whol >> 1;
89: PetscFunctionBeginHot;
90: if (!B && !D && !H) PetscFunctionReturnVoid();
91: for (i = 0; i < nPoints; i++) {
92: const PetscScalar *point = &points[dim * i];
93: PetscReal deta[3] = {0.};
94: PetscReal etaB[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
95: PetscReal etaD[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
96: PetscReal work[8];
98: for (j = 0; j < dim; j++) {
99: PetscReal e, d;
101: e = (PetscRealPart(point[j]) - coordRange[j][0]) / coordRange[j][1];
102: deta[j] = d = 1. / coordRange[j][1];
103: for (k = 0; k < whol; k++) work[k] = etaB[k];
104: for (k = 0; k < half; k++) {
105: etaB[k] = work[2 * k] * e;
106: etaB[k + half] = work[2 * k + 1];
107: }
108: if (H) {
109: for (k = 0; k < whol; k++) work[k] = etaD[k];
110: for (k = 0; k < half; k++) {
111: etaD[k + half] = work[2 * k];
112: etaD[k] = work[2 * k + 1] * d;
113: }
114: }
115: }
116: if (B) {
117: if (datatype == PETSC_SCALAR) {
118: PetscScalar *out = &((PetscScalar *)B)[nc * i];
120: MEdot(out, cf, etaB, (1 << dim), nc, (PetscScalar));
121: } else {
122: PetscReal *out = &((PetscReal *)B)[nc * i];
124: MEdot(out, cf, etaB, (1 << dim), nc, PetscRealPart);
125: }
126: }
127: if (D) {
128: if (datatype == PETSC_SCALAR) {
129: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
131: for (m = 0; m < nc * dim; m++) out[m] = 0.;
132: } else {
133: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
135: for (m = 0; m < nc * dim; m++) out[m] = 0.;
136: }
137: for (j = 0; j < dim; j++) {
138: PetscReal d = deta[j];
140: for (k = 0; k < whol * nc; k++) cfWork[k] = cf[k];
141: for (k = 0; k < whol; k++) work[k] = etaB[k];
142: for (k = 0; k < half; k++) {
143: PetscReal e;
145: etaB[k] = work[2 * k];
146: etaB[k + half] = e = work[2 * k + 1];
148: for (l = 0; l < nc; l++) {
149: cf[k * nc + l] = cfWork[2 * k * nc + l];
150: cf[(k + half) * nc + l] = cfWork[(2 * k + 1) * nc + l];
151: }
152: if (datatype == PETSC_SCALAR) {
153: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
155: for (l = 0; l < nc; l++) out[l * dim + j] += d * e * cf[k * nc + l];
156: } else {
157: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
159: for (l = 0; l < nc; l++) out[l * dim + j] += d * e * PetscRealPart(cf[k * nc + l]);
160: }
161: }
162: }
163: }
164: if (H) {
165: if (datatype == PETSC_SCALAR) {
166: PetscScalar *out = &((PetscScalar *)H)[nc * dim * dim * i];
168: MEHess(out, cf, etaB, etaD, dim, nc, (PetscScalar));
169: } else {
170: PetscReal *out = &((PetscReal *)H)[nc * dim * dim * i];
172: MEHess(out, cf, etaB, etaD, dim, nc, PetscRealPart);
173: }
174: }
175: }
176: PetscFunctionReturnVoid();
177: }
179: static PetscErrorCode DMFieldEvaluate_DA(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
180: {
181: DM dm;
182: DMField_DA *dafield;
183: PetscInt dim;
184: PetscInt N, n, nc;
185: const PetscScalar *array;
186: PetscReal (*coordRange)[2];
188: PetscFunctionBegin;
189: dm = field->dm;
190: nc = field->numComponents;
191: dafield = (DMField_DA *)field->data;
192: PetscCall(DMGetDimension(dm, &dim));
193: PetscCall(VecGetLocalSize(points, &N));
194: PetscCheck(N % dim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Point vector size %" PetscInt_FMT " not divisible by coordinate dimension %" PetscInt_FMT, N, dim);
195: n = N / dim;
196: coordRange = &dafield->coordRange[0];
197: PetscCall(VecGetArrayRead(points, &array));
198: MultilinearEvaluate(dim, coordRange, nc, dafield->cornerCoeffs, dafield->work, n, array, datatype, B, D, H);
199: PetscCall(VecRestoreArrayRead(points, &array));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode DMFieldEvaluateFE_DA(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
204: {
205: PetscInt c, i, j, k, dim, cellsPer[3] = {0}, first[3] = {0}, whol, half;
206: PetscReal stepPer[3] = {0.};
207: PetscReal cellCoordRange[3][2] = {
208: {0., 1.},
209: {0., 1.},
210: {0., 1.}
211: };
212: PetscScalar *cellCoeffs, *work;
213: DM dm;
214: DMDALocalInfo info;
215: PetscInt cStart, cEnd;
216: PetscInt nq, nc;
217: const PetscReal *q;
218: #if PetscDefined(USE_COMPLEX)
219: PetscScalar *qs;
220: #else
221: const PetscScalar *qs;
222: #endif
223: DMField_DA *dafield;
224: PetscBool isStride;
225: const PetscInt *cells = NULL;
226: PetscInt sfirst = -1, stride = -1, nCells;
228: PetscFunctionBegin;
229: dafield = (DMField_DA *)field->data;
230: dm = field->dm;
231: nc = field->numComponents;
232: PetscCall(DMDAGetLocalInfo(dm, &info));
233: dim = info.dim;
234: work = dafield->work;
235: stepPer[0] = 1. / info.mx;
236: stepPer[1] = 1. / info.my;
237: stepPer[2] = 1. / info.mz;
238: first[0] = info.gxs;
239: first[1] = info.gys;
240: first[2] = info.gzs;
241: cellsPer[0] = info.gxm;
242: cellsPer[1] = info.gym;
243: cellsPer[2] = info.gzm;
244: /* TODO: probably take components into account */
245: PetscCall(PetscQuadratureGetData(points, NULL, NULL, &nq, &q, NULL));
246: #if PetscDefined(USE_COMPLEX)
247: PetscCall(DMGetWorkArray(dm, nq * dim, MPIU_SCALAR, &qs));
248: for (i = 0; i < nq * dim; i++) qs[i] = q[i];
249: #else
250: qs = q;
251: #endif
252: PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd));
253: PetscCall(DMGetWorkArray(dm, (1 << dim) * nc, MPIU_SCALAR, &cellCoeffs));
254: whol = (1 << dim);
255: half = whol >> 1;
256: PetscCall(ISGetLocalSize(cellIS, &nCells));
257: PetscCall(PetscObjectTypeCompare((PetscObject)cellIS, ISSTRIDE, &isStride));
258: if (isStride) PetscCall(ISStrideGetInfo(cellIS, &sfirst, &stride));
259: else PetscCall(ISGetIndices(cellIS, &cells));
260: for (c = 0; c < nCells; c++) {
261: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
262: PetscInt rem = cell;
263: PetscInt ijk[3] = {0};
264: void *cB, *cD, *cH;
266: if (datatype == PETSC_SCALAR) {
267: cB = PetscSafePointerPlusOffset((PetscScalar *)B, nc * nq * c);
268: cD = PetscSafePointerPlusOffset((PetscScalar *)D, nc * nq * dim * c);
269: cH = PetscSafePointerPlusOffset((PetscScalar *)H, nc * nq * dim * dim * c);
270: } else {
271: cB = PetscSafePointerPlusOffset((PetscReal *)B, nc * nq * c);
272: cD = PetscSafePointerPlusOffset((PetscReal *)D, nc * nq * dim * c);
273: cH = PetscSafePointerPlusOffset((PetscReal *)H, nc * nq * dim * dim * c);
274: }
275: PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not a cell [%" PetscInt_FMT ",%" PetscInt_FMT "), not implemented yet", cell, cStart, cEnd);
276: for (i = 0; i < nc * whol; i++) work[i] = dafield->cornerCoeffs[i];
277: for (j = 0; j < dim; j++) {
278: PetscReal e, d;
279: ijk[j] = (rem % cellsPer[j]);
280: rem /= cellsPer[j];
282: e = 2. * (ijk[j] + first[j] + 0.5) * stepPer[j] - 1.;
283: d = stepPer[j];
284: for (i = 0; i < half; i++) {
285: for (k = 0; k < nc; k++) {
286: cellCoeffs[i * nc + k] = work[2 * i * nc + k] * d;
287: cellCoeffs[(i + half) * nc + k] = work[2 * i * nc + k] * e + work[(2 * i + 1) * nc + k];
288: }
289: }
290: for (i = 0; i < whol * nc; i++) work[i] = cellCoeffs[i];
291: }
292: MultilinearEvaluate(dim, cellCoordRange, nc, cellCoeffs, dafield->work, nq, qs, datatype, cB, cD, cH);
293: }
294: if (!isStride) PetscCall(ISRestoreIndices(cellIS, &cells));
295: PetscCall(DMRestoreWorkArray(dm, (1 << dim) * nc, MPIU_SCALAR, &cellCoeffs));
296: #if PetscDefined(USE_COMPLEX)
297: PetscCall(DMRestoreWorkArray(dm, nq * dim, MPIU_SCALAR, &qs));
298: #endif
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: static PetscErrorCode DMFieldEvaluateFV_DA(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
303: {
304: PetscInt c, i, dim, cellsPer[3] = {0}, first[3] = {0};
305: PetscReal stepPer[3] = {0.};
306: DM dm;
307: DMDALocalInfo info;
308: PetscInt cStart, cEnd, numCells;
309: PetscInt nc;
310: PetscScalar *points;
311: DMField_DA *dafield;
312: PetscBool isStride;
313: const PetscInt *cells = NULL;
314: PetscInt sfirst = -1, stride = -1;
316: PetscFunctionBegin;
317: dafield = (DMField_DA *)field->data;
318: dm = field->dm;
319: nc = field->numComponents;
320: PetscCall(DMDAGetLocalInfo(dm, &info));
321: dim = info.dim;
322: stepPer[0] = 1. / info.mx;
323: stepPer[1] = 1. / info.my;
324: stepPer[2] = 1. / info.mz;
325: first[0] = info.gxs;
326: first[1] = info.gys;
327: first[2] = info.gzs;
328: cellsPer[0] = info.gxm;
329: cellsPer[1] = info.gym;
330: cellsPer[2] = info.gzm;
331: PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd));
332: PetscCall(ISGetLocalSize(cellIS, &numCells));
333: PetscCall(DMGetWorkArray(dm, dim * numCells, MPIU_SCALAR, &points));
334: PetscCall(PetscObjectTypeCompare((PetscObject)cellIS, ISSTRIDE, &isStride));
335: if (isStride) PetscCall(ISStrideGetInfo(cellIS, &sfirst, &stride));
336: else PetscCall(ISGetIndices(cellIS, &cells));
337: for (c = 0; c < numCells; c++) {
338: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
339: PetscInt rem = cell;
340: PetscInt ijk[3] = {0};
342: PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not a cell [%" PetscInt_FMT ",%" PetscInt_FMT "), not implemented yet", cell, cStart, cEnd);
343: for (i = 0; i < dim; i++) {
344: ijk[i] = (rem % cellsPer[i]);
345: rem /= cellsPer[i];
346: points[dim * c + i] = (ijk[i] + first[i] + 0.5) * stepPer[i];
347: }
348: }
349: if (!isStride) PetscCall(ISRestoreIndices(cellIS, &cells));
350: MultilinearEvaluate(dim, dafield->coordRange, nc, dafield->cornerCoeffs, dafield->work, numCells, points, datatype, B, D, H);
351: PetscCall(DMRestoreWorkArray(dm, dim * numCells, MPIU_SCALAR, &points));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode DMFieldGetDegree_DA(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
356: {
357: DM dm;
358: PetscInt dim, h, imin;
360: PetscFunctionBegin;
361: dm = field->dm;
362: PetscCall(ISGetMinMax(pointIS, &imin, NULL));
363: PetscCall(DMGetDimension(dm, &dim));
364: for (h = 0; h <= dim; h++) {
365: PetscInt hEnd;
367: PetscCall(DMDAGetHeightStratum(dm, h, NULL, &hEnd));
368: if (imin < hEnd) break;
369: }
370: dim -= h;
371: if (minDegree) *minDegree = 1;
372: if (maxDegree) *maxDegree = dim;
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad)
377: {
378: PetscInt h, dim, imax, imin;
379: DM dm;
381: PetscFunctionBegin;
382: dm = field->dm;
383: PetscCall(ISGetMinMax(cellIS, &imin, &imax));
384: PetscCall(DMGetDimension(dm, &dim));
385: *quad = NULL;
386: for (h = 0; h <= dim; h++) {
387: PetscInt hStart, hEnd;
389: PetscCall(DMDAGetHeightStratum(dm, h, &hStart, &hEnd));
390: if (imin >= hStart && imax < hEnd) break;
391: }
392: dim -= h;
393: if (dim > 0) PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static PetscErrorCode DMFieldInitialize_DA(DMField field)
398: {
399: DM dm;
400: Vec coords = NULL;
401: PetscInt dim, i, j, k;
402: DMField_DA *dafield = (DMField_DA *)field->data;
404: PetscFunctionBegin;
405: field->ops->destroy = DMFieldDestroy_DA;
406: field->ops->evaluate = DMFieldEvaluate_DA;
407: field->ops->evaluateFE = DMFieldEvaluateFE_DA;
408: field->ops->evaluateFV = DMFieldEvaluateFV_DA;
409: field->ops->getDegree = DMFieldGetDegree_DA;
410: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA;
411: field->ops->view = DMFieldView_DA;
412: dm = field->dm;
413: PetscCall(DMGetDimension(dm, &dim));
414: PetscCall(DMGetCoordinates(dm, &coords));
415: if (coords) {
416: PetscInt n;
417: const PetscScalar *array;
419: for (j = 0; j < dim; j++) dafield->coordRange[j][0] = dafield->coordRange[j][1] = PETSC_MAX_REAL;
420: PetscCall(VecGetLocalSize(coords, &n));
421: n /= dim;
422: PetscCall(VecGetArrayRead(coords, &array));
423: for (i = 0, k = 0; i < n; i++) {
424: for (j = 0; j < dim; j++, k++) {
425: PetscReal val = PetscRealPart(array[k]);
427: dafield->coordRange[j][0] = PetscMin(dafield->coordRange[j][0], val);
428: dafield->coordRange[j][1] = PetscMin(dafield->coordRange[j][1], -val);
429: }
430: }
431: PetscCall(VecRestoreArrayRead(coords, &array));
432: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &dafield->coordRange[0][0], 2 * dim, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
433: for (j = 0; j < dim; j++) dafield->coordRange[j][1] = -dafield->coordRange[j][1];
434: } else {
435: for (j = 0; j < dim; j++) {
436: dafield->coordRange[j][0] = 0.;
437: dafield->coordRange[j][1] = 1.;
438: }
439: }
440: for (j = 0; j < dim; j++) {
441: PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]);
442: PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]);
444: dafield->coordRange[j][0] = avg;
445: dafield->coordRange[j][1] = dif;
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field)
451: {
452: DMField_DA *dafield;
454: PetscFunctionBegin;
455: PetscCall(PetscNew(&dafield));
456: field->data = dafield;
457: PetscCall(DMFieldInitialize_DA(field));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /*@
462: DMFieldCreateDA - Create a `DMField` of type `DMFIELDDA` that represents a multilinear field on a `DMDA` given by its values at the corners of the reference element.
464: Collective
466: Input Parameters:
467: + dm - the `DMDA` on which the field lives
468: . nc - the number of components of the field
469: - cornerValues - array of length `nc * (1 << dim)` holding the field values at each corner of the reference element, ordered by lexicographic corner index
471: Output Parameter:
472: . field - the newly created `DMField`
474: Level: intermediate
476: Note:
477: Internally the corner values are converted to coefficients of a tensor-product multilinear polynomial so that evaluation at arbitrary points is inexpensive.
479: .seealso: `DMField`, `DMFIELDDA`, `DMDA`, `DMFieldCreate()`, `DMFieldCreateDS()`, `DMFieldCreateShell()`
480: @*/
481: PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues, DMField *field)
482: {
483: DMField b;
484: DMField_DA *dafield;
485: PetscInt dim, nv, i, j, k;
486: PetscInt half;
487: PetscScalar *cv, *cf, *work;
489: PetscFunctionBegin;
490: PetscCall(DMFieldCreate(dm, nc, DMFIELD_VERTEX, &b));
491: PetscCall(DMFieldSetType(b, DMFIELDDA));
492: dafield = (DMField_DA *)b->data;
493: PetscCall(DMGetDimension(dm, &dim));
494: nv = (1 << dim) * nc;
495: PetscCall(PetscMalloc3(nv, &cv, nv, &cf, nv, &work));
496: for (i = 0; i < nv; i++) cv[i] = cornerValues[i];
497: for (i = 0; i < nv; i++) cf[i] = cv[i];
498: dafield->cornerVals = cv;
499: dafield->cornerCoeffs = cf;
500: dafield->work = work;
501: half = (1 << (dim - 1));
502: for (i = 0; i < dim; i++) {
503: PetscScalar *w;
505: w = work;
506: for (j = 0; j < half; j++) {
507: for (k = 0; k < nc; k++) w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]);
508: }
509: w = &work[j * nc];
510: for (j = 0; j < half; j++) {
511: for (k = 0; k < nc; k++) w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]);
512: }
513: for (j = 0; j < nv; j++) cf[j] = work[j];
514: }
515: *field = b;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }