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: {
7: PetscScalar *cornerVals;
8: PetscScalar *cornerCoeffs;
9: PetscScalar *work;
10: PetscReal coordRange[3][2];
11: }
12: DMField_DA;
14: static PetscErrorCode DMFieldDestroy_DA(DMField field)
15: {
16: DMField_DA *dafield;
18: dafield = (DMField_DA *) field->data;
19: PetscFree3(dafield->cornerVals,dafield->cornerCoeffs,dafield->work);
20: PetscFree(dafield);
21: return 0;
22: }
24: static PetscErrorCode DMFieldView_DA(DMField field,PetscViewer viewer)
25: {
26: DMField_DA *dafield = (DMField_DA *) field->data;
27: PetscBool iascii;
29: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
30: if (iascii) {
31: PetscInt i, c, dim;
32: PetscInt nc;
33: DM dm = field->dm;
35: PetscViewerASCIIPrintf(viewer, "Field corner values:\n");
36: PetscViewerASCIIPushTab(viewer);
37: DMGetDimension(dm,&dim);
38: nc = field->numComponents;
39: for (i = 0, c = 0; i < (1 << dim); i++) {
40: PetscInt j;
42: for (j = 0; j < nc; j++, c++) {
43: PetscScalar val = dafield->cornerVals[nc * i + j];
45: #if !defined(PETSC_USE_COMPLEX)
46: PetscViewerASCIIPrintf(viewer,"%g ",(double) val);
47: #else
48: PetscViewerASCIIPrintf(viewer,"%g+i%g ",(double) PetscRealPart(val),(double) PetscImaginaryPart(val));
49: #endif
50: }
51: PetscViewerASCIIPrintf(viewer,"\n");
52: }
53: PetscViewerASCIIPopTab(viewer);
54: }
55: return 0;
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++) { \
64: (y)[_k] += cast((A)[(c) * _l + _k]) * (x)[_l]; \
65: } \
66: } \
67: } while (0)
69: #define MEHess(out,cf,etaB,etaD,dim,nc,cast) \
70: do { \
71: PetscInt _m, _j, _k; \
72: for (_m = 0; _m < (nc) * (dim) * (dim); _m++) (out)[_m] = 0.; \
73: for (_j = 0; _j < (dim); _j++) { \
74: for (_k = _j + 1; _k < (dim); _k++) { \
75: PetscInt _ind = (1 << _j) + (1 << _k); \
76: for (_m = 0; _m < (nc); _m++) { \
77: PetscScalar c = (cf)[_m] * (etaB)[_ind] * (etaD)[_ind]; \
78: (out)[(_m * (dim) + _k) * (dim) + _j] += cast(c); \
79: (out)[(_m * (dim) + _j) * (dim) + _k] += cast(c); \
80: } \
81: } \
82: } \
83: } while (0)
85: 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)
86: {
87: PetscInt i, j, k, l, m;
88: PetscInt whol = 1 << dim;
89: PetscInt half = whol >> 1;
92: if (!B && !D && !H) return;
93: for (i = 0; i < nPoints; i++) {
94: const PetscScalar *point = &points[dim * i];
95: PetscReal deta[3] = {0.};
96: PetscReal etaB[8] = {1.,1.,1.,1.,1.,1.,1.,1.};
97: PetscReal etaD[8] = {1.,1.,1.,1.,1.,1.,1.,1.};
98: PetscReal work[8];
100: for (j = 0; j < dim; j++) {
101: PetscReal e, d;
103: e = (PetscRealPart(point[j]) - coordRange[j][0]) / coordRange[j][1];
104: deta[j] = d = 1. / coordRange[j][1];
105: for (k = 0; k < whol; k++) {work[k] = etaB[k];}
106: for (k = 0; k < half; k++) {
107: etaB[k] = work[2 * k] * e;
108: etaB[k + half] = work[2 * k + 1];
109: }
110: if (H) {
111: for (k = 0; k < whol; k++) {work[k] = etaD[k];}
112: for (k = 0; k < half; k++) {
113: etaD[k + half] = work[2 * k];
114: etaD[k ] = work[2 * k + 1] * d;
115: }
116: }
117: }
118: if (B) {
119: if (datatype == PETSC_SCALAR) {
120: PetscScalar *out = &((PetscScalar *)B)[nc * i];
122: MEdot(out,cf,etaB,(1 << dim),nc,(PetscScalar));
123: } else {
124: PetscReal *out = &((PetscReal *)B)[nc * i];
126: MEdot(out,cf,etaB,(1 << dim),nc,PetscRealPart);
127: }
128: }
129: if (D) {
130: if (datatype == PETSC_SCALAR) {
131: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
133: for (m = 0; m < nc * dim; m++) out[m] = 0.;
134: } else {
135: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
137: for (m = 0; m < nc * dim; m++) out[m] = 0.;
138: }
139: for (j = 0; j < dim; j++) {
140: PetscReal d = deta[j];
142: for (k = 0; k < whol * nc; k++) {cfWork[k] = cf[k];}
143: for (k = 0; k < whol; k++) {work[k] = etaB[k];}
144: for (k = 0; k < half; k++) {
145: PetscReal e;
147: etaB[k] = work[2 * k];
148: etaB[k + half] = e = work[2 * k + 1];
150: for (l = 0; l < nc; l++) {
151: cf[ k * nc + l] = cfWork[ 2 * k * nc + l];
152: cf[(k + half) * nc + l] = cfWork[(2 * k + 1) * nc + l];
153: }
154: if (datatype == PETSC_SCALAR) {
155: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
157: for (l = 0; l < nc; l++) {
158: out[l * dim + j] += d * e * cf[k * nc + l];
159: }
160: } else {
161: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
163: for (l = 0; l < nc; l++) {
164: out[l * dim + j] += d * e * PetscRealPart(cf[k * nc + l]);
165: }
166: }
167: }
168: }
169: }
170: if (H) {
171: if (datatype == PETSC_SCALAR) {
172: PetscScalar *out = &((PetscScalar *)H)[nc * dim * dim * i];
174: MEHess(out,cf,etaB,etaD,dim,nc,(PetscScalar));
175: } else {
176: PetscReal *out = &((PetscReal *)H)[nc * dim * dim * i];
178: MEHess(out,cf,etaB,etaD,dim,nc,PetscRealPart);
179: }
180: }
181: }
182: return;
183: }
185: static PetscErrorCode DMFieldEvaluate_DA(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
186: {
187: DM dm;
188: DMField_DA *dafield;
189: PetscInt dim;
190: PetscInt N, n, nc;
191: const PetscScalar *array;
192: PetscReal (*coordRange)[2];
194: dm = field->dm;
195: nc = field->numComponents;
196: dafield = (DMField_DA *) field->data;
197: DMGetDimension(dm,&dim);
198: VecGetLocalSize(points,&N);
200: n = N / dim;
201: coordRange = &(dafield->coordRange[0]);
202: VecGetArrayRead(points,&array);
203: MultilinearEvaluate(dim,coordRange,nc,dafield->cornerCoeffs,dafield->work,n,array,datatype,B,D,H);
204: VecRestoreArrayRead(points,&array);
205: return 0;
206: }
208: static PetscErrorCode DMFieldEvaluateFE_DA(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
209: {
210: PetscInt c, i, j, k, dim, cellsPer[3] = {0}, first[3] = {0}, whol, half;
211: PetscReal stepPer[3] = {0.};
212: PetscReal cellCoordRange[3][2] = {{0.,1.},{0.,1.},{0.,1.}};
213: PetscScalar *cellCoeffs, *work;
214: DM dm;
215: DMDALocalInfo info;
216: PetscInt cStart, cEnd;
217: PetscInt nq, nc;
218: const PetscReal *q;
219: #if defined(PETSC_USE_COMPLEX)
220: PetscScalar *qs;
221: #else
222: const PetscScalar *qs;
223: #endif
224: DMField_DA *dafield;
225: PetscBool isStride;
226: const PetscInt *cells = NULL;
227: PetscInt sfirst = -1, stride = -1, nCells;
229: dafield = (DMField_DA *) field->data;
230: dm = field->dm;
231: nc = field->numComponents;
232: 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: PetscQuadratureGetData(points, NULL, NULL, &nq, &q, NULL);
246: #if defined(PETSC_USE_COMPLEX)
247: 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: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
253: DMGetWorkArray(dm,(1 << dim) * nc,MPIU_SCALAR,&cellCoeffs);
254: whol = (1 << dim);
255: half = whol >> 1;
256: ISGetLocalSize(cellIS,&nCells);
257: PetscObjectTypeCompare((PetscObject)cellIS,ISSTRIDE,&isStride);
258: if (isStride) {
259: ISStrideGetInfo(cellIS,&sfirst,&stride);
260: } else {
261: ISGetIndices(cellIS,&cells);
262: }
263: for (c = 0; c < nCells; c++) {
264: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
265: PetscInt rem = cell;
266: PetscInt ijk[3] = {0};
267: void *cB, *cD, *cH;
269: if (datatype == PETSC_SCALAR) {
270: cB = B ? &((PetscScalar *)B)[nc * nq * c] : NULL;
271: cD = D ? &((PetscScalar *)D)[nc * nq * dim * c] : NULL;
272: cH = H ? &((PetscScalar *)H)[nc * nq * dim * dim * c] : NULL;
273: } else {
274: cB = B ? &((PetscReal *)B)[nc * nq * c] : NULL;
275: cD = D ? &((PetscReal *)D)[nc * nq * dim * c] : NULL;
276: cH = H ? &((PetscReal *)H)[nc * nq * dim * dim * c] : NULL;
277: }
279: for (i = 0; i < nc * whol; i++) {work[i] = dafield->cornerCoeffs[i];}
280: for (j = 0; j < dim; j++) {
281: PetscReal e, d;
282: ijk[j] = (rem % cellsPer[j]);
283: rem /= cellsPer[j];
285: e = 2. * (ijk[j] + first[j] + 0.5) * stepPer[j] - 1.;
286: d = stepPer[j];
287: for (i = 0; i < half; i++) {
288: for (k = 0; k < nc; k++) {
289: cellCoeffs[ i * nc + k] = work[ 2 * i * nc + k] * d;
290: cellCoeffs[(i + half) * nc + k] = work[ 2 * i * nc + k] * e + work[(2 * i + 1) * nc + k];
291: }
292: }
293: for (i = 0; i < whol * nc; i++) {work[i] = cellCoeffs[i];}
294: }
295: MultilinearEvaluate(dim,cellCoordRange,nc,cellCoeffs,dafield->work,nq,qs,datatype,cB,cD,cH);
296: }
297: if (!isStride) {
298: ISRestoreIndices(cellIS,&cells);
299: }
300: DMRestoreWorkArray(dm,(1 << dim) * nc,MPIU_SCALAR,&cellCoeffs);
301: #if defined(PETSC_USE_COMPLEX)
302: DMRestoreWorkArray(dm,nq * dim,MPIU_SCALAR,&qs);
303: #endif
304: return 0;
305: }
307: static PetscErrorCode DMFieldEvaluateFV_DA(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
308: {
309: PetscInt c, i, dim, cellsPer[3] = {0}, first[3] = {0};
310: PetscReal stepPer[3] = {0.};
311: DM dm;
312: DMDALocalInfo info;
313: PetscInt cStart, cEnd, numCells;
314: PetscInt nc;
315: PetscScalar *points;
316: DMField_DA *dafield;
317: PetscBool isStride;
318: const PetscInt *cells = NULL;
319: PetscInt sfirst = -1, stride = -1;
321: dafield = (DMField_DA *) field->data;
322: dm = field->dm;
323: nc = field->numComponents;
324: DMDAGetLocalInfo(dm,&info);
325: dim = info.dim;
326: stepPer[0] = 1./ info.mx;
327: stepPer[1] = 1./ info.my;
328: stepPer[2] = 1./ info.mz;
329: first[0] = info.gxs;
330: first[1] = info.gys;
331: first[2] = info.gzs;
332: cellsPer[0] = info.gxm;
333: cellsPer[1] = info.gym;
334: cellsPer[2] = info.gzm;
335: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
336: ISGetLocalSize(cellIS,&numCells);
337: DMGetWorkArray(dm,dim * numCells,MPIU_SCALAR,&points);
338: PetscObjectTypeCompare((PetscObject)cellIS,ISSTRIDE,&isStride);
339: if (isStride) {
340: ISStrideGetInfo(cellIS,&sfirst,&stride);
341: } else {
342: ISGetIndices(cellIS,&cells);
343: }
344: for (c = 0; c < numCells; c++) {
345: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
346: PetscInt rem = cell;
347: PetscInt ijk[3] = {0};
350: for (i = 0; i < dim; i++) {
351: ijk[i] = (rem % cellsPer[i]);
352: rem /= cellsPer[i];
353: points[dim * c + i] = (ijk[i] + first[i] + 0.5) * stepPer[i];
354: }
355: }
356: if (!isStride) {
357: ISRestoreIndices(cellIS,&cells);
358: }
359: MultilinearEvaluate(dim,dafield->coordRange,nc,dafield->cornerCoeffs,dafield->work,numCells,points,datatype,B,D,H);
360: DMRestoreWorkArray(dm,dim * numCells,MPIU_SCALAR,&points);
361: return 0;
362: }
364: static PetscErrorCode DMFieldGetDegree_DA(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
365: {
366: DM dm;
367: PetscInt dim, h, imin;
369: dm = field->dm;
370: ISGetMinMax(pointIS,&imin,NULL);
371: DMGetDimension(dm,&dim);
372: for (h = 0; h <= dim; h++) {
373: PetscInt hEnd;
375: DMDAGetHeightStratum(dm,h,NULL,&hEnd);
376: if (imin < hEnd) break;
377: }
378: dim -= h;
379: if (minDegree) *minDegree = 1;
380: if (maxDegree) *maxDegree = dim;
381: return 0;
382: }
384: static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad)
385: {
386: PetscInt h, dim, imax, imin;
387: DM dm;
389: dm = field->dm;
390: ISGetMinMax(cellIS,&imin,&imax);
391: DMGetDimension(dm,&dim);
392: *quad = NULL;
393: for (h = 0; h <= dim; h++) {
394: PetscInt hStart, hEnd;
396: DMDAGetHeightStratum(dm,h,&hStart,&hEnd);
397: if (imin >= hStart && imax < hEnd) break;
398: }
399: dim -= h;
400: if (dim > 0) {
401: PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);
402: }
404: return 0;
405: }
407: static PetscErrorCode DMFieldInitialize_DA(DMField field)
408: {
409: DM dm;
410: Vec coords = NULL;
411: PetscInt dim, i, j, k;
412: DMField_DA *dafield = (DMField_DA *) field->data;
414: field->ops->destroy = DMFieldDestroy_DA;
415: field->ops->evaluate = DMFieldEvaluate_DA;
416: field->ops->evaluateFE = DMFieldEvaluateFE_DA;
417: field->ops->evaluateFV = DMFieldEvaluateFV_DA;
418: field->ops->getDegree = DMFieldGetDegree_DA;
419: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA;
420: field->ops->view = DMFieldView_DA;
421: dm = field->dm;
422: DMGetDimension(dm,&dim);
423: if (dm->coordinates) coords = dm->coordinates;
424: else if (dm->coordinatesLocal) coords = dm->coordinatesLocal;
425: if (coords) {
426: PetscInt n;
427: const PetscScalar *array;
428: PetscReal mins[3][2] = {{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL}};
430: VecGetLocalSize(coords,&n);
431: n /= dim;
432: VecGetArrayRead(coords,&array);
433: for (i = 0, k = 0; i < n; i++) {
434: for (j = 0; j < dim; j++, k++) {
435: PetscReal val = PetscRealPart(array[k]);
437: mins[j][0] = PetscMin(mins[j][0],val);
438: mins[j][1] = PetscMin(mins[j][1],-val);
439: }
440: }
441: VecRestoreArrayRead(coords,&array);
442: MPIU_Allreduce((PetscReal *) mins,&(dafield->coordRange[0][0]),2*dim,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)dm));
443: for (j = 0; j < dim; j++) {
444: dafield->coordRange[j][1] = -dafield->coordRange[j][1];
445: }
446: } else {
447: for (j = 0; j < dim; j++) {
448: dafield->coordRange[j][0] = 0.;
449: dafield->coordRange[j][1] = 1.;
450: }
451: }
452: for (j = 0; j < dim; j++) {
453: PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]);
454: PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]);
456: dafield->coordRange[j][0] = avg;
457: dafield->coordRange[j][1] = dif;
458: }
459: return 0;
460: }
462: PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field)
463: {
464: DMField_DA *dafield;
466: PetscNewLog(field,&dafield);
467: field->data = dafield;
468: DMFieldInitialize_DA(field);
469: return 0;
470: }
472: PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues,DMField *field)
473: {
474: DMField b;
475: DMField_DA *dafield;
476: PetscInt dim, nv, i, j, k;
477: PetscInt half;
478: PetscScalar *cv, *cf, *work;
480: DMFieldCreate(dm,nc,DMFIELD_VERTEX,&b);
481: DMFieldSetType(b,DMFIELDDA);
482: dafield = (DMField_DA *) b->data;
483: DMGetDimension(dm,&dim);
484: nv = (1 << dim) * nc;
485: PetscMalloc3(nv,&cv,nv,&cf,nv,&work);
486: for (i = 0; i < nv; i++) cv[i] = cornerValues[i];
487: for (i = 0; i < nv; i++) cf[i] = cv[i];
488: dafield->cornerVals = cv;
489: dafield->cornerCoeffs = cf;
490: dafield->work = work;
491: half = (1 << (dim - 1));
492: for (i = 0; i < dim; i++) {
493: PetscScalar *w;
495: w = work;
496: for (j = 0; j < half; j++) {
497: for (k = 0; k < nc; k++) {
498: w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]);
499: }
500: }
501: w = &work[j * nc];
502: for (j = 0; j < half; j++) {
503: for (k = 0; k < nc; k++) {
504: w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]);
505: }
506: }
507: for (j = 0; j < nv; j++) {cf[j] = work[j];}
508: }
509: *field = b;
510: return 0;
511: }