Actual source code: ex2.c

  1: static char help[] = "Interpolation Tests for Plex\n\n";

  3: #include <petscsnes.h>
  4: #include <petscdmplex.h>
  5: #include <petscdmda.h>
  6: #include <petscds.h>

  8: typedef enum {
  9:   CENTROID,
 10:   GRID,
 11:   GRID_REPLICATED
 12: } PointType;

 14: typedef struct {
 15:   PointType pointType; /* Point generation mechanism */
 16: } AppCtx;

 18: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 19: {
 20:   PetscInt d, c;

 23:   for (c = 0; c < Nc; ++c) {
 24:     u[c] = 0.0;
 25:     for (d = 0; d < dim; ++d) u[c] += x[d];
 26:   }
 27:   return 0;
 28: }

 30: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 31: {
 32:   const char *pointTypes[3] = {"centroid", "grid", "grid_replicated"};
 33:   PetscInt    pt;

 35:   options->pointType = CENTROID;
 36:   PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");
 37:   pt = options->pointType;
 38:   PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL);
 39:   options->pointType = (PointType)pt;
 40:   PetscOptionsEnd();
 41:   return 0;
 42: }

 44: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
 45: {
 46:   DMCreate(comm, dm);
 47:   DMSetType(*dm, DMPLEX);
 48:   DMSetFromOptions(*dm);
 49:   DMViewFromOptions(*dm, NULL, "-dm_view");
 50:   return 0;
 51: }

 53: static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
 54: {
 55:   PetscSection coordSection;
 56:   Vec          coordsLocal;
 57:   PetscInt     spaceDim, p;
 58:   PetscMPIInt  rank;

 60:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 61:   DMGetCoordinatesLocal(dm, &coordsLocal);
 62:   DMGetCoordinateSection(dm, &coordSection);
 63:   DMGetCoordinateDim(dm, &spaceDim);
 64:   DMPlexGetHeightStratum(dm, 0, NULL, Np);
 65:   PetscCalloc1(*Np * spaceDim, pcoords);
 66:   for (p = 0; p < *Np; ++p) {
 67:     PetscScalar *coords = NULL;
 68:     PetscInt     size, num, n, d;

 70:     DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords);
 71:     num = size / spaceDim;
 72:     for (n = 0; n < num; ++n) {
 73:       for (d = 0; d < spaceDim; ++d) (*pcoords)[p * spaceDim + d] += PetscRealPart(coords[n * spaceDim + d]) / num;
 74:     }
 75:     PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, p);
 76:     for (d = 0; d < spaceDim; ++d) {
 77:       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p * spaceDim + d]);
 78:       if (d < spaceDim - 1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");
 79:     }
 80:     PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");
 81:     DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords);
 82:   }
 83:   PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
 84:   *pointsAllProcs = PETSC_FALSE;
 85:   return 0;
 86: }

 88: static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
 89: {
 90:   DM            da;
 91:   DMDALocalInfo info;
 92:   PetscInt      N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
 93:   PetscReal    *h;
 94:   PetscMPIInt   rank;

 96:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 97:   DMGetDimension(dm, &dim);
 98:   DMGetCoordinateDim(dm, &spaceDim);
 99:   PetscCalloc1(spaceDim, &ind);
100:   PetscCalloc1(spaceDim, &h);
101:   h[0] = 1.0 / (N - 1);
102:   h[1] = 1.0 / (N - 1);
103:   h[2] = 1.0 / (N - 1);
104:   DMDACreate(PetscObjectComm((PetscObject)dm), &da);
105:   DMSetDimension(da, dim);
106:   DMDASetSizes(da, N, N, N);
107:   DMDASetDof(da, 1);
108:   DMDASetStencilWidth(da, 1);
109:   DMSetUp(da);
110:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
111:   DMDAGetLocalInfo(da, &info);
112:   *Np = info.xm * info.ym * info.zm;
113:   PetscCalloc1(*Np * spaceDim, pcoords);
114:   for (k = info.zs; k < info.zs + info.zm; ++k) {
115:     ind[2] = k;
116:     for (j = info.ys; j < info.ys + info.ym; ++j) {
117:       ind[1] = j;
118:       for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
119:         ind[0] = i;

121:         for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
122:         PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n);
123:         for (d = 0; d < spaceDim; ++d) {
124:           PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]);
125:           if (d < spaceDim - 1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");
126:         }
127:         PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");
128:       }
129:     }
130:   }
131:   DMDestroy(&da);
132:   PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
133:   PetscFree(ind);
134:   PetscFree(h);
135:   *pointsAllProcs = PETSC_FALSE;
136:   return 0;
137: }

139: static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
140: {
141:   PetscInt    N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
142:   PetscReal  *h;
143:   PetscMPIInt rank;

146:   DMGetDimension(dm, &dim);
147:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
148:   DMGetCoordinateDim(dm, &spaceDim);
149:   PetscCalloc1(spaceDim, &ind);
150:   PetscCalloc1(spaceDim, &h);
151:   h[0] = 1.0 / (N - 1);
152:   h[1] = 1.0 / (N - 1);
153:   h[2] = 1.0 / (N - 1);
154:   *Np  = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1);
155:   PetscCalloc1(*Np * spaceDim, pcoords);
156:   for (k = 0; k < N; ++k) {
157:     ind[2] = k;
158:     for (j = 0; j < N; ++j) {
159:       ind[1] = j;
160:       for (i = 0; i < N; ++i, ++n) {
161:         ind[0] = i;

163:         for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
164:         PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n);
165:         for (d = 0; d < spaceDim; ++d) {
166:           PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]);
167:           if (d < spaceDim - 1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");
168:         }
169:         PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");
170:       }
171:     }
172:   }
173:   PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
174:   *pointsAllProcs = PETSC_TRUE;
175:   PetscFree(ind);
176:   PetscFree(h);
177:   return 0;
178: }

180: static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
181: {
182:   *pointsAllProcs = PETSC_FALSE;
183:   switch (ctx->pointType) {
184:   case CENTROID:
185:     CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx);
186:     break;
187:   case GRID:
188:     CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx);
189:     break;
190:   case GRID_REPLICATED:
191:     CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx);
192:     break;
193:   default:
194:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int)ctx->pointType);
195:   }
196:   return 0;
197: }

199: int main(int argc, char **argv)
200: {
201:   AppCtx ctx;
202:   PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
203:   DM                  dm;
204:   PetscFE             fe;
205:   DMInterpolationInfo interpolator;
206:   Vec                 lu, fieldVals;
207:   PetscScalar        *vals;
208:   const PetscScalar  *ivals, *vcoords;
209:   PetscReal          *pcoords;
210:   PetscBool           simplex, pointsAllProcs = PETSC_TRUE;
211:   PetscInt            dim, spaceDim, Nc, c, Np, p;
212:   PetscMPIInt         rank, size;
213:   PetscViewer         selfviewer;

216:   PetscInitialize(&argc, &argv, NULL, help);
217:   ProcessOptions(PETSC_COMM_WORLD, &ctx);
218:   CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);
219:   DMGetDimension(dm, &dim);
220:   DMGetCoordinateDim(dm, &spaceDim);
221:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
222:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
223:   /* Create points */
224:   CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx);
225:   /* Create interpolator */
226:   DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator);
227:   DMInterpolationSetDim(interpolator, spaceDim);
228:   DMInterpolationAddPoints(interpolator, Np, pcoords);
229:   DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE);
230:   /* Check locations */
231:   for (c = 0; c < interpolator->n; ++c) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " is in Cell %" PetscInt_FMT "\n", rank, c, interpolator->cells[c]);
232:   PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
233:   VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD);
234:   /* Setup Discretization */
235:   Nc = dim;
236:   DMPlexIsSimplex(dm, &simplex);
237:   PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, Nc, simplex, NULL, -1, &fe);
238:   DMSetField(dm, 0, NULL, (PetscObject)fe);
239:   DMCreateDS(dm);
240:   PetscFEDestroy(&fe);
241:   /* Create function */
242:   PetscCalloc2(Nc, &funcs, Nc, &vals);
243:   for (c = 0; c < Nc; ++c) funcs[c] = linear;
244:   DMGetLocalVector(dm, &lu);
245:   DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu);
246:   PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
247:   PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer);
248:   PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank);
249:   VecView(lu, selfviewer);
250:   PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer);
251:   PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
252:   PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
253:   /* Check interpolant */
254:   VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals);
255:   DMInterpolationSetDof(interpolator, Nc);
256:   DMInterpolationEvaluate(interpolator, dm, lu, fieldVals);
257:   for (p = 0; p < size; ++p) {
258:     if (p == rank) {
259:       PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank);
260:       VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF);
261:     }
262:     PetscBarrier((PetscObject)dm);
263:   }
264:   VecGetArrayRead(interpolator->coords, &vcoords);
265:   VecGetArrayRead(fieldVals, &ivals);
266:   for (p = 0; p < interpolator->n; ++p) {
267:     for (c = 0; c < Nc; ++c) {
268: #if defined(PETSC_USE_COMPLEX)
269:       PetscReal vcoordsReal[3];
270:       PetscInt  i;

272:       for (i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
273: #else
274:       const PetscReal *vcoordsReal = &vcoords[p * spaceDim];
275: #endif
276:       (*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL);
277:       if (PetscAbsScalar(ivals[p * Nc + c] - vals[c]) > PETSC_SQRT_MACHINE_EPSILON)
278:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid interpolated value %g != %g (%" PetscInt_FMT ", %" PetscInt_FMT ")", (double)PetscRealPart(ivals[p * Nc + c]), (double)PetscRealPart(vals[c]), p, c);
279:     }
280:   }
281:   VecRestoreArrayRead(interpolator->coords, &vcoords);
282:   VecRestoreArrayRead(fieldVals, &ivals);
283:   /* Cleanup */
284:   PetscFree(pcoords);
285:   PetscFree2(funcs, vals);
286:   VecDestroy(&fieldVals);
287:   DMRestoreLocalVector(dm, &lu);
288:   DMInterpolationDestroy(&interpolator);
289:   DMDestroy(&dm);
290:   PetscFinalize();
291:   return 0;
292: }

294: /*TEST

296:   testset:
297:     requires: ctetgen
298:     args: -dm_plex_dim 3 -petscspace_degree 1

300:     test:
301:       suffix: 0
302:     test:
303:       suffix: 1
304:       args: -dm_refine 2
305:     test:
306:       suffix: 2
307:       nsize: 2
308:       args: -petscpartitioner_type simple
309:     test:
310:       suffix: 3
311:       nsize: 2
312:       args: -dm_refine 2 -petscpartitioner_type simple
313:     test:
314:       suffix: 4
315:       nsize: 5
316:       args: -petscpartitioner_type simple
317:     test:
318:       suffix: 5
319:       nsize: 5
320:       args: -dm_refine 2 -petscpartitioner_type simple
321:     test:
322:       suffix: 6
323:       args: -point_type grid
324:     test:
325:       suffix: 7
326:       args: -dm_refine 2 -point_type grid
327:     test:
328:       suffix: 8
329:       nsize: 2
330:       args: -petscpartitioner_type simple -point_type grid
331:     test:
332:       suffix: 9
333:       args: -point_type grid_replicated
334:     test:
335:       suffix: 10
336:       nsize: 2
337:       args: -petscpartitioner_type simple -point_type grid_replicated
338:     test:
339:       suffix: 11
340:       nsize: 2
341:       args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated

343: TEST*/