Actual source code: dmfieldda.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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: }