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