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