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: }