Actual source code: dmfieldda.c
petsc-3.13.6 2020-09-29
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;
20: dafield = (DMField_DA *) field->data;
21: PetscFree3(dafield->cornerVals,dafield->cornerCoeffs,dafield->work);
22: PetscFree(dafield);
23: return(0);
24: }
26: static PetscErrorCode DMFieldView_DA(DMField field,PetscViewer viewer)
27: {
28: DMField_DA *dafield = (DMField_DA *) field->data;
29: PetscBool iascii;
33: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
34: if (iascii) {
35: PetscInt i, c, dim;
36: PetscInt nc;
37: DM dm = field->dm;
39: PetscViewerASCIIPrintf(viewer, "Field corner values:\n");
40: PetscViewerASCIIPushTab(viewer);
41: DMGetDimension(dm,&dim);
42: nc = field->numComponents;
43: for (i = 0, c = 0; i < (1 << dim); i++) {
44: PetscInt j;
46: for (j = 0; j < nc; j++, c++) {
47: PetscScalar val = dafield->cornerVals[nc * i + j];
49: #if !defined(PETSC_USE_COMPLEX)
50: PetscViewerASCIIPrintf(viewer,"%g ",(double) val);
51: #else
52: PetscViewerASCIIPrintf(viewer,"%g+i%g ",(double) PetscRealPart(val),(double) PetscImaginaryPart(val));
53: #endif
54: }
55: PetscViewerASCIIPrintf(viewer,"\n");
56: }
57: PetscViewerASCIIPopTab(viewer);
58: }
59: return(0);
60: }
62: #define MEdot(y,A,x,m,c,cast) \
63: do { \
64: PetscInt _k, _l; \
65: for (_k = 0; _k < (c); _k++) (y)[_k] = 0.; \
66: for (_l = 0; _l < (m); _l++) { \
67: for (_k = 0; _k < (c); _k++) { \
68: (y)[_k] += cast((A)[(c) * _l + _k]) * (x)[_l]; \
69: } \
70: } \
71: } while (0)
73: #define MEHess(out,cf,etaB,etaD,dim,nc,cast) \
74: do { \
75: PetscInt _m, _j, _k; \
76: for (_m = 0; _m < (nc) * (dim) * (dim); _m++) (out)[_m] = 0.; \
77: for (_j = 0; _j < (dim); _j++) { \
78: for (_k = _j + 1; _k < (dim); _k++) { \
79: PetscInt _ind = (1 << _j) + (1 << _k); \
80: for (_m = 0; _m < (nc); _m++) { \
81: PetscScalar c = (cf)[_m] * (etaB)[_ind] * (etaD)[_ind]; \
82: (out)[(_m * (dim) + _k) * (dim) + _j] += cast(c); \
83: (out)[(_m * (dim) + _j) * (dim) + _k] += cast(c); \
84: } \
85: } \
86: } \
87: } while (0)
89: 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)
90: {
91: PetscInt i, j, k, l, m;
92: PetscInt whol = 1 << dim;
93: PetscInt half = whol >> 1;
96: if (!B && !D && !H) PetscFunctionReturnVoid();
97: for (i = 0; i < nPoints; i++) {
98: const PetscScalar *point = &points[dim * i];
99: PetscReal deta[3] = {0.};
100: PetscReal etaB[8] = {1.,1.,1.,1.,1.,1.,1.,1.};
101: PetscReal etaD[8] = {1.,1.,1.,1.,1.,1.,1.,1.};
102: PetscReal work[8];
104: for (j = 0; j < dim; j++) {
105: PetscReal e, d;
107: e = (PetscRealPart(point[j]) - coordRange[j][0]) / coordRange[j][1];
108: deta[j] = d = 1. / coordRange[j][1];
109: for (k = 0; k < whol; k++) {work[k] = etaB[k];}
110: for (k = 0; k < half; k++) {
111: etaB[k] = work[2 * k] * e;
112: etaB[k + half] = work[2 * k + 1];
113: }
114: if (H) {
115: for (k = 0; k < whol; k++) {work[k] = etaD[k];}
116: for (k = 0; k < half; k++) {
117: etaD[k + half] = work[2 * k];
118: etaD[k ] = work[2 * k + 1] * d;
119: }
120: }
121: }
122: if (B) {
123: if (datatype == PETSC_SCALAR) {
124: PetscScalar *out = &((PetscScalar *)B)[nc * i];
126: MEdot(out,cf,etaB,(1 << dim),nc,(PetscScalar));
127: } else {
128: PetscReal *out = &((PetscReal *)B)[nc * i];
130: MEdot(out,cf,etaB,(1 << dim),nc,PetscRealPart);
131: }
132: }
133: if (D) {
134: if (datatype == PETSC_SCALAR) {
135: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
137: for (m = 0; m < nc * dim; m++) out[m] = 0.;
138: } else {
139: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
141: for (m = 0; m < nc * dim; m++) out[m] = 0.;
142: }
143: for (j = 0; j < dim; j++) {
144: PetscReal d = deta[j];
146: for (k = 0; k < whol * nc; k++) {cfWork[k] = cf[k];}
147: for (k = 0; k < whol; k++) {work[k] = etaB[k];}
148: for (k = 0; k < half; k++) {
149: PetscReal e;
151: etaB[k] = work[2 * k];
152: etaB[k + half] = e = work[2 * k + 1];
154: for (l = 0; l < nc; l++) {
155: cf[ k * nc + l] = cfWork[ 2 * k * nc + l];
156: cf[(k + half) * nc + l] = cfWork[(2 * k + 1) * nc + l];
157: }
158: if (datatype == PETSC_SCALAR) {
159: PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
161: for (l = 0; l < nc; l++) {
162: out[l * dim + j] += d * e * cf[k * nc + l];
163: }
164: } else {
165: PetscReal *out = &((PetscReal *)D)[nc * dim * i];
167: for (l = 0; l < nc; l++) {
168: out[l * dim + j] += d * e * PetscRealPart(cf[k * nc + l]);
169: }
170: }
171: }
172: }
173: }
174: if (H) {
175: if (datatype == PETSC_SCALAR) {
176: PetscScalar *out = &((PetscScalar *)H)[nc * dim * dim * i];
178: MEHess(out,cf,etaB,etaD,dim,nc,(PetscScalar));
179: } else {
180: PetscReal *out = &((PetscReal *)H)[nc * dim * dim * i];
182: MEHess(out,cf,etaB,etaD,dim,nc,PetscRealPart);
183: }
184: }
185: }
186: PetscFunctionReturnVoid();
187: }
189: static PetscErrorCode DMFieldEvaluate_DA(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
190: {
191: DM dm;
192: DMField_DA *dafield;
193: PetscInt dim;
194: PetscInt N, n, nc;
195: const PetscScalar *array;
196: PetscReal (*coordRange)[2];
200: dm = field->dm;
201: nc = field->numComponents;
202: dafield = (DMField_DA *) field->data;
203: DMGetDimension(dm,&dim);
204: VecGetLocalSize(points,&N);
205: if (N % dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point vector size %D not divisible by coordinate dimension %D\n",N,dim);
206: n = N / dim;
207: coordRange = &(dafield->coordRange[0]);
208: VecGetArrayRead(points,&array);
209: MultilinearEvaluate(dim,coordRange,nc,dafield->cornerCoeffs,dafield->work,n,array,datatype,B,D,H);
210: VecRestoreArrayRead(points,&array);
211: return(0);
212: }
214: static PetscErrorCode DMFieldEvaluateFE_DA(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
215: {
216: PetscInt c, i, j, k, dim, cellsPer[3] = {0}, first[3] = {0}, whol, half;
217: PetscReal stepPer[3] = {0.};
218: PetscReal cellCoordRange[3][2] = {{0.,1.},{0.,1.},{0.,1.}};
219: PetscScalar *cellCoeffs, *work;
220: DM dm;
221: DMDALocalInfo info;
222: PetscInt cStart, cEnd;
223: PetscInt nq, nc;
224: const PetscReal *q;
225: #if defined(PETSC_USE_COMPLEX)
226: PetscScalar *qs;
227: #else
228: const PetscScalar *qs;
229: #endif
230: DMField_DA *dafield;
231: PetscBool isStride;
232: const PetscInt *cells = NULL;
233: PetscInt sfirst = -1, stride = -1, nCells;
237: dafield = (DMField_DA *) field->data;
238: dm = field->dm;
239: nc = field->numComponents;
240: DMDAGetLocalInfo(dm,&info);
241: dim = info.dim;
242: work = dafield->work;
243: stepPer[0] = 1./ info.mx;
244: stepPer[1] = 1./ info.my;
245: stepPer[2] = 1./ info.mz;
246: first[0] = info.gxs;
247: first[1] = info.gys;
248: first[2] = info.gzs;
249: cellsPer[0] = info.gxm;
250: cellsPer[1] = info.gym;
251: cellsPer[2] = info.gzm;
252: /* TODO: probably take components into account */
253: PetscQuadratureGetData(points, NULL, NULL, &nq, &q, NULL);
254: #if defined(PETSC_USE_COMPLEX)
255: DMGetWorkArray(dm,nq * dim,MPIU_SCALAR,&qs);
256: for (i = 0; i < nq * dim; i++) qs[i] = q[i];
257: #else
258: qs = q;
259: #endif
260: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
261: DMGetWorkArray(dm,(1 << dim) * nc,MPIU_SCALAR,&cellCoeffs);
262: whol = (1 << dim);
263: half = whol >> 1;
264: ISGetLocalSize(cellIS,&nCells);
265: PetscObjectTypeCompare((PetscObject)cellIS,ISSTRIDE,&isStride);
266: if (isStride) {
267: ISStrideGetInfo(cellIS,&sfirst,&stride);
268: } else {
269: ISGetIndices(cellIS,&cells);
270: }
271: for (c = 0; c < nCells; c++) {
272: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
273: PetscInt rem = cell;
274: PetscInt ijk[3] = {0};
275: void *cB, *cD, *cH;
277: if (datatype == PETSC_SCALAR) {
278: cB = B ? &((PetscScalar *)B)[nc * nq * c] : NULL;
279: cD = D ? &((PetscScalar *)D)[nc * nq * dim * c] : NULL;
280: cH = H ? &((PetscScalar *)H)[nc * nq * dim * dim * c] : NULL;
281: } else {
282: cB = B ? &((PetscReal *)B)[nc * nq * c] : NULL;
283: cD = D ? &((PetscReal *)D)[nc * nq * dim * c] : NULL;
284: cH = H ? &((PetscReal *)H)[nc * nq * dim * dim * c] : NULL;
285: }
286: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Point %D not a cell [%D,%D), not implemented yet",cell,cStart,cEnd);
287: for (i = 0; i < nc * whol; i++) {work[i] = dafield->cornerCoeffs[i];}
288: for (j = 0; j < dim; j++) {
289: PetscReal e, d;
290: ijk[j] = (rem % cellsPer[j]);
291: rem /= cellsPer[j];
293: e = 2. * (ijk[j] + first[j] + 0.5) * stepPer[j] - 1.;
294: d = stepPer[j];
295: for (i = 0; i < half; i++) {
296: for (k = 0; k < nc; k++) {
297: cellCoeffs[ i * nc + k] = work[ 2 * i * nc + k] * d;
298: cellCoeffs[(i + half) * nc + k] = work[ 2 * i * nc + k] * e + work[(2 * i + 1) * nc + k];
299: }
300: }
301: for (i = 0; i < whol * nc; i++) {work[i] = cellCoeffs[i];}
302: }
303: MultilinearEvaluate(dim,cellCoordRange,nc,cellCoeffs,dafield->work,nq,qs,datatype,cB,cD,cH);
304: }
305: if (!isStride) {
306: ISRestoreIndices(cellIS,&cells);
307: }
308: DMRestoreWorkArray(dm,(1 << dim) * nc,MPIU_SCALAR,&cellCoeffs);
309: #if defined(PETSC_USE_COMPLEX)
310: DMRestoreWorkArray(dm,nq * dim,MPIU_SCALAR,&qs);
311: #endif
312: return(0);
313: }
315: static PetscErrorCode DMFieldEvaluateFV_DA(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
316: {
317: PetscInt c, i, dim, cellsPer[3] = {0}, first[3] = {0};
318: PetscReal stepPer[3] = {0.};
319: DM dm;
320: DMDALocalInfo info;
321: PetscInt cStart, cEnd, numCells;
322: PetscInt nc;
323: PetscScalar *points;
324: DMField_DA *dafield;
325: PetscBool isStride;
326: const PetscInt *cells = NULL;
327: PetscInt sfirst = -1, stride = -1;
331: dafield = (DMField_DA *) field->data;
332: dm = field->dm;
333: nc = field->numComponents;
334: DMDAGetLocalInfo(dm,&info);
335: dim = info.dim;
336: stepPer[0] = 1./ info.mx;
337: stepPer[1] = 1./ info.my;
338: stepPer[2] = 1./ info.mz;
339: first[0] = info.gxs;
340: first[1] = info.gys;
341: first[2] = info.gzs;
342: cellsPer[0] = info.gxm;
343: cellsPer[1] = info.gym;
344: cellsPer[2] = info.gzm;
345: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
346: ISGetLocalSize(cellIS,&numCells);
347: DMGetWorkArray(dm,dim * numCells,MPIU_SCALAR,&points);
348: PetscObjectTypeCompare((PetscObject)cellIS,ISSTRIDE,&isStride);
349: if (isStride) {
350: ISStrideGetInfo(cellIS,&sfirst,&stride);
351: } else {
352: ISGetIndices(cellIS,&cells);
353: }
354: for (c = 0; c < numCells; c++) {
355: PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
356: PetscInt rem = cell;
357: PetscInt ijk[3] = {0};
359: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Point %D not a cell [%D,%D), not implemented yet",cell,cStart,cEnd);
360: for (i = 0; i < dim; i++) {
361: ijk[i] = (rem % cellsPer[i]);
362: rem /= cellsPer[i];
363: points[dim * c + i] = (ijk[i] + first[i] + 0.5) * stepPer[i];
364: }
365: }
366: if (!isStride) {
367: ISRestoreIndices(cellIS,&cells);
368: }
369: MultilinearEvaluate(dim,dafield->coordRange,nc,dafield->cornerCoeffs,dafield->work,numCells,points,datatype,B,D,H);
370: DMRestoreWorkArray(dm,dim * numCells,MPIU_SCALAR,&points);
371: return(0);
372: }
374: static PetscErrorCode DMFieldGetDegree_DA(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
375: {
376: DM dm;
377: PetscInt dim, h, imin;
381: dm = field->dm;
382: ISGetMinMax(pointIS,&imin,NULL);
383: DMGetDimension(dm,&dim);
384: for (h = 0; h <= dim; h++) {
385: PetscInt hEnd;
387: DMDAGetHeightStratum(dm,h,NULL,&hEnd);
388: if (imin < hEnd) break;
389: }
390: dim -= h;
391: if (minDegree) *minDegree = 1;
392: if (maxDegree) *maxDegree = dim;
393: return(0);
394: }
396: static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad)
397: {
398: PetscInt h, dim, imax, imin;
399: DM dm;
403: dm = field->dm;
404: ISGetMinMax(cellIS,&imin,&imax);
405: DMGetDimension(dm,&dim);
406: *quad = NULL;
407: for (h = 0; h <= dim; h++) {
408: PetscInt hStart, hEnd;
410: DMDAGetHeightStratum(dm,h,&hStart,&hEnd);
411: if (imin >= hStart && imax < hEnd) break;
412: }
413: dim -= h;
414: if (dim > 0) {
415: PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);
416: }
418: return(0);
419: }
421: static PetscErrorCode DMFieldInitialize_DA(DMField field)
422: {
423: DM dm;
424: Vec coords = NULL;
425: PetscInt dim, i, j, k;
426: DMField_DA *dafield = (DMField_DA *) field->data;
430: field->ops->destroy = DMFieldDestroy_DA;
431: field->ops->evaluate = DMFieldEvaluate_DA;
432: field->ops->evaluateFE = DMFieldEvaluateFE_DA;
433: field->ops->evaluateFV = DMFieldEvaluateFV_DA;
434: field->ops->getDegree = DMFieldGetDegree_DA;
435: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA;
436: field->ops->view = DMFieldView_DA;
437: dm = field->dm;
438: DMGetDimension(dm,&dim);
439: if (dm->coordinates) coords = dm->coordinates;
440: else if (dm->coordinatesLocal) coords = dm->coordinatesLocal;
441: if (coords) {
442: PetscInt n;
443: const PetscScalar *array;
444: PetscReal mins[3][2] = {{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL}};
446: VecGetLocalSize(coords,&n);
447: n /= dim;
448: VecGetArrayRead(coords,&array);
449: for (i = 0, k = 0; i < n; i++) {
450: for (j = 0; j < dim; j++, k++) {
451: PetscReal val = PetscRealPart(array[k]);
453: mins[j][0] = PetscMin(mins[j][0],val);
454: mins[j][1] = PetscMin(mins[j][1],-val);
455: }
456: }
457: VecRestoreArrayRead(coords,&array);
458: MPIU_Allreduce((PetscReal *) mins,&(dafield->coordRange[0][0]),2*dim,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)dm));
459: for (j = 0; j < dim; j++) {
460: dafield->coordRange[j][1] = -dafield->coordRange[j][1];
461: }
462: } else {
463: for (j = 0; j < dim; j++) {
464: dafield->coordRange[j][0] = 0.;
465: dafield->coordRange[j][1] = 1.;
466: }
467: }
468: for (j = 0; j < dim; j++) {
469: PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]);
470: PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]);
472: dafield->coordRange[j][0] = avg;
473: dafield->coordRange[j][1] = dif;
474: }
475: return(0);
476: }
478: PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field)
479: {
480: DMField_DA *dafield;
484: PetscNewLog(field,&dafield);
485: field->data = dafield;
486: DMFieldInitialize_DA(field);
487: return(0);
488: }
490: PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues,DMField *field)
491: {
492: DMField b;
493: DMField_DA *dafield;
494: PetscInt dim, nv, i, j, k;
495: PetscInt half;
496: PetscScalar *cv, *cf, *work;
500: DMFieldCreate(dm,nc,DMFIELD_VERTEX,&b);
501: DMFieldSetType(b,DMFIELDDA);
502: dafield = (DMField_DA *) b->data;
503: DMGetDimension(dm,&dim);
504: nv = (1 << dim) * nc;
505: PetscMalloc3(nv,&cv,nv,&cf,nv,&work);
506: for (i = 0; i < nv; i++) cv[i] = cornerValues[i];
507: for (i = 0; i < nv; i++) cf[i] = cv[i];
508: dafield->cornerVals = cv;
509: dafield->cornerCoeffs = cf;
510: dafield->work = work;
511: half = (1 << (dim - 1));
512: for (i = 0; i < dim; i++) {
513: PetscScalar *w;
515: w = work;
516: for (j = 0; j < half; j++) {
517: for (k = 0; k < nc; k++) {
518: w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]);
519: }
520: }
521: w = &work[j * nc];
522: for (j = 0; j < half; j++) {
523: for (k = 0; k < nc; k++) {
524: w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]);
525: }
526: }
527: for (j = 0; j < nv; j++) {cf[j] = work[j];}
528: }
529: *field = b;
530: return(0);
531: }