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