1 static char help[] = "Interpolation Tests for Plex\n\n";
2
3 #include <petscsnes.h>
4 #include <petscdmplex.h>
5 #include <petscdmda.h>
6 #include <petscds.h>
7
8 typedef enum {
9 CENTROID,
10 GRID,
11 GRID_REPLICATED
12 } PointType;
13
14 typedef enum {
15 CONSTANT,
16 LINEAR
17 } FuncType;
18
19 typedef struct {
20 PointType pointType; // Point generation mechanism
21 FuncType funcType; // Type of interpolated function
22 PetscBool useFV; // Use finite volume, instead of finite element
23 } AppCtx;
24
constant(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)25 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
26 {
27 PetscFunctionBeginUser;
28 for (PetscInt c = 0; c < Nc; ++c) u[c] = c + 1.;
29 PetscFunctionReturn(PETSC_SUCCESS);
30 }
31
linear(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)32 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
33 {
34 PetscFunctionBeginUser;
35 PetscCheck(Nc == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Something is wrong: %" PetscInt_FMT, Nc);
36 for (PetscInt c = 0; c < Nc; ++c) {
37 u[c] = 0.0;
38 for (PetscInt d = 0; d < dim; ++d) u[c] += x[d];
39 }
40 PetscFunctionReturn(PETSC_SUCCESS);
41 }
42
ProcessOptions(MPI_Comm comm,AppCtx * options)43 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
44 {
45 const char *pointTypes[3] = {"centroid", "grid", "grid_replicated"};
46 const char *funcTypes[2] = {"constant", "linear"};
47 PetscInt pt, fn;
48
49 PetscFunctionBegin;
50 options->pointType = CENTROID;
51 options->funcType = LINEAR;
52 options->useFV = PETSC_FALSE;
53
54 PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");
55 pt = options->pointType;
56 PetscCall(PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL));
57 options->pointType = (PointType)pt;
58 fn = options->funcType;
59 PetscCall(PetscOptionsEList("-func_type", "The function type", "ex2.c", funcTypes, 2, funcTypes[options->funcType], &fn, NULL));
60 options->funcType = (FuncType)fn;
61 PetscCall(PetscOptionsBool("-use_fv", "Use finite volumes, instead of finite elements", "ex2.c", options->useFV, &options->useFV, NULL));
62 PetscOptionsEnd();
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
CreateMesh(MPI_Comm comm,AppCtx * ctx,DM * dm)66 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
67 {
68 PetscFunctionBegin;
69 PetscCall(DMCreate(comm, dm));
70 PetscCall(DMSetType(*dm, DMPLEX));
71 PetscCall(DMSetFromOptions(*dm));
72 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
CreatePoints_Centroid(DM dm,PetscInt * Np,PetscReal ** pcoords,PetscBool * pointsAllProcs,AppCtx * ctx)76 static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
77 {
78 PetscSection coordSection;
79 Vec coordsLocal;
80 PetscInt spaceDim, p;
81 PetscMPIInt rank;
82
83 PetscFunctionBegin;
84 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
85 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
86 PetscCall(DMGetCoordinateSection(dm, &coordSection));
87 PetscCall(DMGetCoordinateDim(dm, &spaceDim));
88 PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, Np));
89 PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
90 for (p = 0; p < *Np; ++p) {
91 PetscScalar *coords = NULL;
92 PetscInt size, num, n, d;
93
94 PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords));
95 num = size / spaceDim;
96 for (n = 0; n < num; ++n) {
97 for (d = 0; d < spaceDim; ++d) (*pcoords)[p * spaceDim + d] += PetscRealPart(coords[n * spaceDim + d]) / num;
98 }
99 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, p));
100 for (d = 0; d < spaceDim; ++d) {
101 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p * spaceDim + d]));
102 if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
103 }
104 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
105 PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords));
106 }
107 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
108 *pointsAllProcs = PETSC_FALSE;
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
CreatePoints_Grid(DM dm,PetscInt * Np,PetscReal ** pcoords,PetscBool * pointsAllProcs,AppCtx * ctx)112 static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
113 {
114 DM da;
115 DMDALocalInfo info;
116 PetscInt N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
117 PetscReal *h;
118 PetscMPIInt rank;
119
120 PetscFunctionBegin;
121 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
122 PetscCall(DMGetDimension(dm, &dim));
123 PetscCall(DMGetCoordinateDim(dm, &spaceDim));
124 PetscCall(PetscCalloc1(spaceDim, &ind));
125 PetscCall(PetscCalloc1(spaceDim, &h));
126 h[0] = 1.0 / (N - 1);
127 h[1] = 1.0 / (N - 1);
128 h[2] = 1.0 / (N - 1);
129 PetscCall(DMDACreate(PetscObjectComm((PetscObject)dm), &da));
130 PetscCall(DMSetDimension(da, dim));
131 PetscCall(DMDASetSizes(da, N, N, N));
132 PetscCall(DMDASetDof(da, 1));
133 PetscCall(DMDASetStencilWidth(da, 1));
134 PetscCall(DMSetUp(da));
135 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
136 PetscCall(DMDAGetLocalInfo(da, &info));
137 *Np = info.xm * info.ym * info.zm;
138 PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
139 for (k = info.zs; k < info.zs + info.zm; ++k) {
140 ind[2] = k;
141 for (j = info.ys; j < info.ys + info.ym; ++j) {
142 ind[1] = j;
143 for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
144 ind[0] = i;
145
146 for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
147 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
148 for (d = 0; d < spaceDim; ++d) {
149 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
150 if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
151 }
152 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
153 }
154 }
155 }
156 PetscCall(DMDestroy(&da));
157 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
158 PetscCall(PetscFree(ind));
159 PetscCall(PetscFree(h));
160 *pointsAllProcs = PETSC_FALSE;
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
CreatePoints_GridReplicated(DM dm,PetscInt * Np,PetscReal ** pcoords,PetscBool * pointsAllProcs,AppCtx * ctx)164 static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
165 {
166 PetscInt N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
167 PetscReal *h;
168 PetscMPIInt rank;
169
170 PetscFunctionBeginUser;
171 PetscCall(DMGetDimension(dm, &dim));
172 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
173 PetscCall(DMGetCoordinateDim(dm, &spaceDim));
174 PetscCall(PetscCalloc1(spaceDim, &ind));
175 PetscCall(PetscCalloc1(spaceDim, &h));
176 h[0] = 1.0 / (N - 1);
177 h[1] = 1.0 / (N - 1);
178 h[2] = 1.0 / (N - 1);
179 *Np = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1);
180 PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
181 for (k = 0; k < N; ++k) {
182 ind[2] = k;
183 for (j = 0; j < N; ++j) {
184 ind[1] = j;
185 for (i = 0; i < N; ++i, ++n) {
186 ind[0] = i;
187
188 for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
189 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
190 for (d = 0; d < spaceDim; ++d) {
191 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
192 if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
193 }
194 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
195 }
196 }
197 }
198 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
199 *pointsAllProcs = PETSC_TRUE;
200 PetscCall(PetscFree(ind));
201 PetscCall(PetscFree(h));
202 PetscFunctionReturn(PETSC_SUCCESS);
203 }
204
CreatePoints(DM dm,PetscInt * Np,PetscReal ** pcoords,PetscBool * pointsAllProcs,AppCtx * ctx)205 static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
206 {
207 PetscFunctionBegin;
208 *pointsAllProcs = PETSC_FALSE;
209 switch (ctx->pointType) {
210 case CENTROID:
211 PetscCall(CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx));
212 break;
213 case GRID:
214 PetscCall(CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx));
215 break;
216 case GRID_REPLICATED:
217 PetscCall(CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx));
218 break;
219 default:
220 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int)ctx->pointType);
221 }
222 PetscFunctionReturn(PETSC_SUCCESS);
223 }
224
CreateDiscretization(DM dm,PetscInt Nc,AppCtx * ctx)225 static PetscErrorCode CreateDiscretization(DM dm, PetscInt Nc, AppCtx *ctx)
226 {
227 PetscFunctionBegin;
228 if (ctx->useFV) {
229 PetscFV fv;
230 PetscInt cdim;
231
232 PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)dm), &fv));
233 PetscCall(PetscObjectSetName((PetscObject)fv, "phi"));
234 PetscCall(PetscFVSetFromOptions(fv));
235 PetscCall(PetscFVSetNumComponents(fv, Nc));
236 PetscCall(DMGetCoordinateDim(dm, &cdim));
237 PetscCall(PetscFVSetSpatialDimension(fv, cdim));
238 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv));
239 PetscCall(PetscFVDestroy(&fv));
240 PetscCall(DMCreateDS(dm));
241 } else {
242 PetscFE fe;
243 DMPolytopeType ct;
244 PetscInt dim, cStart;
245
246 PetscCall(DMGetDimension(dm, &dim));
247 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
248 PetscCall(DMPlexGetCellType(dm, cStart, &ct));
249 PetscCall(PetscFECreateByCell(PetscObjectComm((PetscObject)dm), dim, Nc, ct, NULL, -1, &fe));
250 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
251 PetscCall(PetscFEDestroy(&fe));
252 PetscCall(DMCreateDS(dm));
253 }
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
main(int argc,char ** argv)257 int main(int argc, char **argv)
258 {
259 AppCtx ctx;
260 PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
261 DM dm;
262 DMInterpolationInfo interpolator;
263 Vec lu, fieldVals;
264 PetscScalar *vals;
265 const PetscScalar *ivals, *vcoords;
266 PetscReal *pcoords;
267 PetscBool pointsAllProcs = PETSC_TRUE;
268 PetscInt dim, spaceDim, Nc, c, Np, p;
269 PetscMPIInt rank, size;
270 PetscViewer selfviewer;
271
272 PetscFunctionBeginUser;
273 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
274 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
275 PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
276 PetscCall(DMGetDimension(dm, &dim));
277 PetscCall(DMGetCoordinateDim(dm, &spaceDim));
278 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
279 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
280 /* Create points */
281 PetscCall(CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx));
282 /* Create interpolator */
283 PetscCall(DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator));
284 PetscCall(DMInterpolationSetDim(interpolator, spaceDim));
285 PetscCall(DMInterpolationAddPoints(interpolator, Np, pcoords));
286 PetscCall(DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE));
287 /* Check locations */
288 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]));
289 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
290 PetscCall(VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD));
291 Nc = dim;
292 PetscCall(CreateDiscretization(dm, Nc, &ctx));
293 /* Create function */
294 PetscCall(PetscCalloc2(Nc, &funcs, Nc, &vals));
295 switch (ctx.funcType) {
296 case CONSTANT:
297 for (c = 0; c < Nc; ++c) funcs[c] = constant;
298 break;
299 case LINEAR:
300 for (c = 0; c < Nc; ++c) funcs[c] = linear;
301 break;
302 default:
303 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid function type: %d", (int)ctx.funcType);
304 }
305 PetscCall(DMGetLocalVector(dm, &lu));
306 PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu));
307 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
308 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
309 PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank));
310 PetscCall(VecView(lu, selfviewer));
311 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
312 PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
313 /* Check interpolant */
314 PetscCall(VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals));
315 PetscCall(DMInterpolationSetDof(interpolator, Nc));
316 PetscCall(DMInterpolationEvaluate(interpolator, dm, lu, fieldVals));
317 for (p = 0; p < size; ++p) {
318 if (p == rank) {
319 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank));
320 PetscCall(VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF));
321 }
322 PetscCall(PetscBarrier((PetscObject)dm));
323 }
324 PetscCall(VecGetArrayRead(interpolator->coords, &vcoords));
325 PetscCall(VecGetArrayRead(fieldVals, &ivals));
326 for (p = 0; p < interpolator->n; ++p) {
327 for (c = 0; c < Nc; ++c) {
328 #if defined(PETSC_USE_COMPLEX)
329 PetscReal vcoordsReal[3];
330
331 for (PetscInt i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
332 #else
333 const PetscReal *vcoordsReal = &vcoords[p * spaceDim];
334 #endif
335 PetscCall((*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL));
336 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);
337 }
338 }
339 PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords));
340 PetscCall(VecRestoreArrayRead(fieldVals, &ivals));
341 /* Cleanup */
342 PetscCall(PetscFree(pcoords));
343 PetscCall(PetscFree2(funcs, vals));
344 PetscCall(VecDestroy(&fieldVals));
345 PetscCall(DMRestoreLocalVector(dm, &lu));
346 PetscCall(DMInterpolationDestroy(&interpolator));
347 PetscCall(DMDestroy(&dm));
348 PetscCall(PetscFinalize());
349 return 0;
350 }
351
352 /*TEST
353
354 testset:
355 requires: ctetgen
356 args: -dm_plex_dim 3 -petscspace_degree 1
357
358 test:
359 suffix: 0
360 test:
361 suffix: 1
362 args: -dm_refine 2
363 test:
364 suffix: 2
365 nsize: 2
366 args: -petscpartitioner_type simple
367 test:
368 suffix: 3
369 nsize: 2
370 args: -dm_refine 2 -petscpartitioner_type simple
371 test:
372 suffix: 4
373 nsize: 5
374 args: -petscpartitioner_type simple
375 test:
376 suffix: 5
377 nsize: 5
378 args: -dm_refine 2 -petscpartitioner_type simple
379 test:
380 suffix: 6
381 args: -point_type grid
382 test:
383 suffix: 7
384 args: -dm_refine 2 -point_type grid
385 test:
386 suffix: 8
387 nsize: 2
388 args: -petscpartitioner_type simple -point_type grid
389 test:
390 suffix: 9
391 args: -point_type grid_replicated
392 test:
393 suffix: 10
394 nsize: 2
395 args: -petscpartitioner_type simple -point_type grid_replicated
396 test:
397 suffix: 11
398 nsize: 2
399 args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated
400
401 test:
402 suffix: fv_0
403 requires: triangle
404 args: -use_fv -func_type constant
405
406 TEST*/
407