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*/