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;
22: PetscFunctionBeginUser;
23: PetscCheck(Nc == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Something is wrong: %" PetscInt_FMT, Nc);
24: for (c = 0; c < Nc; ++c) {
25: u[c] = 0.0;
26: for (d = 0; d < dim; ++d) u[c] += x[d];
27: }
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
32: {
33: const char *pointTypes[3] = {"centroid", "grid", "grid_replicated"};
34: PetscInt pt;
36: PetscFunctionBegin;
37: options->pointType = CENTROID;
38: PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");
39: pt = options->pointType;
40: PetscCall(PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL));
41: options->pointType = (PointType)pt;
42: PetscOptionsEnd();
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
47: {
48: PetscFunctionBegin;
49: PetscCall(DMCreate(comm, dm));
50: PetscCall(DMSetType(*dm, DMPLEX));
51: PetscCall(DMSetFromOptions(*dm));
52: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
57: {
58: PetscSection coordSection;
59: Vec coordsLocal;
60: PetscInt spaceDim, p;
61: PetscMPIInt rank;
63: PetscFunctionBegin;
64: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
65: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
66: PetscCall(DMGetCoordinateSection(dm, &coordSection));
67: PetscCall(DMGetCoordinateDim(dm, &spaceDim));
68: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, Np));
69: PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
70: for (p = 0; p < *Np; ++p) {
71: PetscScalar *coords = NULL;
72: PetscInt size, num, n, d;
74: PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords));
75: num = size / spaceDim;
76: for (n = 0; n < num; ++n) {
77: for (d = 0; d < spaceDim; ++d) (*pcoords)[p * spaceDim + d] += PetscRealPart(coords[n * spaceDim + d]) / num;
78: }
79: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, p));
80: for (d = 0; d < spaceDim; ++d) {
81: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p * spaceDim + d]));
82: if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
83: }
84: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
85: PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords));
86: }
87: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
88: *pointsAllProcs = PETSC_FALSE;
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
93: {
94: DM da;
95: DMDALocalInfo info;
96: PetscInt N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
97: PetscReal *h;
98: PetscMPIInt rank;
100: PetscFunctionBegin;
101: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
102: PetscCall(DMGetDimension(dm, &dim));
103: PetscCall(DMGetCoordinateDim(dm, &spaceDim));
104: PetscCall(PetscCalloc1(spaceDim, &ind));
105: PetscCall(PetscCalloc1(spaceDim, &h));
106: h[0] = 1.0 / (N - 1);
107: h[1] = 1.0 / (N - 1);
108: h[2] = 1.0 / (N - 1);
109: PetscCall(DMDACreate(PetscObjectComm((PetscObject)dm), &da));
110: PetscCall(DMSetDimension(da, dim));
111: PetscCall(DMDASetSizes(da, N, N, N));
112: PetscCall(DMDASetDof(da, 1));
113: PetscCall(DMDASetStencilWidth(da, 1));
114: PetscCall(DMSetUp(da));
115: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
116: PetscCall(DMDAGetLocalInfo(da, &info));
117: *Np = info.xm * info.ym * info.zm;
118: PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
119: for (k = info.zs; k < info.zs + info.zm; ++k) {
120: ind[2] = k;
121: for (j = info.ys; j < info.ys + info.ym; ++j) {
122: ind[1] = j;
123: for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
124: ind[0] = i;
126: for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
127: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
128: for (d = 0; d < spaceDim; ++d) {
129: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
130: if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
131: }
132: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
133: }
134: }
135: }
136: PetscCall(DMDestroy(&da));
137: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
138: PetscCall(PetscFree(ind));
139: PetscCall(PetscFree(h));
140: *pointsAllProcs = PETSC_FALSE;
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
145: {
146: PetscInt N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
147: PetscReal *h;
148: PetscMPIInt rank;
150: PetscFunctionBeginUser;
151: PetscCall(DMGetDimension(dm, &dim));
152: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
153: PetscCall(DMGetCoordinateDim(dm, &spaceDim));
154: PetscCall(PetscCalloc1(spaceDim, &ind));
155: PetscCall(PetscCalloc1(spaceDim, &h));
156: h[0] = 1.0 / (N - 1);
157: h[1] = 1.0 / (N - 1);
158: h[2] = 1.0 / (N - 1);
159: *Np = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1);
160: PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
161: for (k = 0; k < N; ++k) {
162: ind[2] = k;
163: for (j = 0; j < N; ++j) {
164: ind[1] = j;
165: for (i = 0; i < N; ++i, ++n) {
166: ind[0] = i;
168: for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
169: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
170: for (d = 0; d < spaceDim; ++d) {
171: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
172: if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
173: }
174: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
175: }
176: }
177: }
178: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
179: *pointsAllProcs = PETSC_TRUE;
180: PetscCall(PetscFree(ind));
181: PetscCall(PetscFree(h));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
186: {
187: PetscFunctionBegin;
188: *pointsAllProcs = PETSC_FALSE;
189: switch (ctx->pointType) {
190: case CENTROID:
191: PetscCall(CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx));
192: break;
193: case GRID:
194: PetscCall(CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx));
195: break;
196: case GRID_REPLICATED:
197: PetscCall(CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx));
198: break;
199: default:
200: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int)ctx->pointType);
201: }
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: int main(int argc, char **argv)
206: {
207: AppCtx ctx;
208: PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
209: DM dm;
210: PetscFE fe;
211: DMInterpolationInfo interpolator;
212: Vec lu, fieldVals;
213: PetscScalar *vals;
214: const PetscScalar *ivals, *vcoords;
215: PetscReal *pcoords;
216: PetscBool simplex, pointsAllProcs = PETSC_TRUE;
217: PetscInt dim, spaceDim, Nc, c, Np, p;
218: PetscMPIInt rank, size;
219: PetscViewer selfviewer;
221: PetscFunctionBeginUser;
222: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
223: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
224: PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
225: PetscCall(DMGetDimension(dm, &dim));
226: PetscCall(DMGetCoordinateDim(dm, &spaceDim));
227: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
228: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
229: /* Create points */
230: PetscCall(CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx));
231: /* Create interpolator */
232: PetscCall(DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator));
233: PetscCall(DMInterpolationSetDim(interpolator, spaceDim));
234: PetscCall(DMInterpolationAddPoints(interpolator, Np, pcoords));
235: PetscCall(DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE));
236: /* Check locations */
237: for (c = 0; c < interpolator->n; ++c) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " is in Cell %" PetscInt_FMT "\n", rank, c, interpolator->cells[c]));
238: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
239: PetscCall(VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD));
240: /* Setup Discretization */
241: Nc = dim;
242: PetscCall(DMPlexIsSimplex(dm, &simplex));
243: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, Nc, simplex, NULL, -1, &fe));
244: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
245: PetscCall(DMCreateDS(dm));
246: PetscCall(PetscFEDestroy(&fe));
247: /* Create function */
248: PetscCall(PetscCalloc2(Nc, &funcs, Nc, &vals));
249: for (c = 0; c < Nc; ++c) funcs[c] = linear;
250: PetscCall(DMGetLocalVector(dm, &lu));
251: PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu));
252: PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
253: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
254: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank));
255: PetscCall(VecView(lu, selfviewer));
256: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
257: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
258: PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
259: /* Check interpolant */
260: PetscCall(VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals));
261: PetscCall(DMInterpolationSetDof(interpolator, Nc));
262: PetscCall(DMInterpolationEvaluate(interpolator, dm, lu, fieldVals));
263: for (p = 0; p < size; ++p) {
264: if (p == rank) {
265: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank));
266: PetscCall(VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF));
267: }
268: PetscCall(PetscBarrier((PetscObject)dm));
269: }
270: PetscCall(VecGetArrayRead(interpolator->coords, &vcoords));
271: PetscCall(VecGetArrayRead(fieldVals, &ivals));
272: for (p = 0; p < interpolator->n; ++p) {
273: for (c = 0; c < Nc; ++c) {
274: #if defined(PETSC_USE_COMPLEX)
275: PetscReal vcoordsReal[3];
277: for (PetscInt i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
278: #else
279: const PetscReal *vcoordsReal = &vcoords[p * spaceDim];
280: #endif
281: PetscCall((*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL));
282: PetscCheck(PetscAbsScalar(ivals[p * Nc + c] - vals[c]) <= PETSC_SQRT_MACHINE_EPSILON, 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);
283: }
284: }
285: PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords));
286: PetscCall(VecRestoreArrayRead(fieldVals, &ivals));
287: /* Cleanup */
288: PetscCall(PetscFree(pcoords));
289: PetscCall(PetscFree2(funcs, vals));
290: PetscCall(VecDestroy(&fieldVals));
291: PetscCall(DMRestoreLocalVector(dm, &lu));
292: PetscCall(DMInterpolationDestroy(&interpolator));
293: PetscCall(DMDestroy(&dm));
294: PetscCall(PetscFinalize());
295: return 0;
296: }
298: /*TEST
300: testset:
301: requires: ctetgen
302: args: -dm_plex_dim 3 -petscspace_degree 1
304: test:
305: suffix: 0
306: test:
307: suffix: 1
308: args: -dm_refine 2
309: test:
310: suffix: 2
311: nsize: 2
312: args: -petscpartitioner_type simple
313: test:
314: suffix: 3
315: nsize: 2
316: args: -dm_refine 2 -petscpartitioner_type simple
317: test:
318: suffix: 4
319: nsize: 5
320: args: -petscpartitioner_type simple
321: test:
322: suffix: 5
323: nsize: 5
324: args: -dm_refine 2 -petscpartitioner_type simple
325: test:
326: suffix: 6
327: args: -point_type grid
328: test:
329: suffix: 7
330: args: -dm_refine 2 -point_type grid
331: test:
332: suffix: 8
333: nsize: 2
334: args: -petscpartitioner_type simple -point_type grid
335: test:
336: suffix: 9
337: args: -point_type grid_replicated
338: test:
339: suffix: 10
340: nsize: 2
341: args: -petscpartitioner_type simple -point_type grid_replicated
342: test:
343: suffix: 11
344: nsize: 2
345: args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated
347: TEST*/