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