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 PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 15 { 16 PetscInt d, c; 17 18 PetscCheckFalse(Nc != 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Something is wrong: %D", Nc); 19 for (c = 0; c < Nc; ++c) { 20 u[c] = 0.0; 21 for (d = 0; d < dim; ++d) u[c] += x[d]; 22 } 23 return 0; 24 } 25 26 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 27 { 28 const char *pointTypes[3] = {"centroid", "grid", "grid_replicated"}; 29 PetscInt pt; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 options->pointType = CENTROID; 34 35 ierr = PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");CHKERRQ(ierr); 36 pt = options->pointType; 37 ierr = PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL);CHKERRQ(ierr); 38 options->pointType = (PointType) pt; 39 ierr = PetscOptionsEnd();CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) 44 { 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 49 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 50 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 51 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx) 56 { 57 PetscSection coordSection; 58 Vec coordsLocal; 59 PetscInt spaceDim, p; 60 PetscMPIInt rank; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 65 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 66 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 67 ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr); 68 ierr = DMPlexGetHeightStratum(dm, 0, NULL, Np);CHKERRQ(ierr); 69 ierr = PetscCalloc1(*Np * spaceDim, pcoords);CHKERRQ(ierr); 70 for (p = 0; p < *Np; ++p) { 71 PetscScalar *coords = NULL; 72 PetscInt size, num, n, d; 73 74 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords);CHKERRQ(ierr); 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 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D (", rank, p);CHKERRQ(ierr); 80 for (d = 0; d < spaceDim; ++d) { 81 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p*spaceDim+d]);CHKERRQ(ierr); 82 if (d < spaceDim-1) {ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");CHKERRQ(ierr);} 83 } 84 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");CHKERRQ(ierr); 85 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords);CHKERRQ(ierr); 86 } 87 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr); 88 *pointsAllProcs = PETSC_FALSE; 89 PetscFunctionReturn(0); 90 } 91 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; 99 PetscErrorCode ierr; 100 101 PetscFunctionBegin; 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: SETERRQ(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, Nc, 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 Nc = dim; 237 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 238 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, NULL, -1, &fe);CHKERRQ(ierr); 239 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 240 ierr = DMCreateDS(dm);CHKERRQ(ierr); 241 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 242 /* Create function */ 243 ierr = PetscCalloc2(Nc, &funcs, Nc, &vals);CHKERRQ(ierr); 244 for (c = 0; c < Nc; ++c) funcs[c] = linear; 245 ierr = DMGetLocalVector(dm, &lu);CHKERRQ(ierr); 246 ierr = DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu);CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 248 ierr = PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&selfviewer);CHKERRQ(ierr); 249 ierr = PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank);CHKERRQ(ierr); 250 ierr = VecView(lu,selfviewer);CHKERRQ(ierr); 251 ierr = PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&selfviewer);CHKERRQ(ierr); 252 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 253 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 254 /* Check interpolant */ 255 ierr = VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals);CHKERRQ(ierr); 256 ierr = DMInterpolationSetDof(interpolator, Nc);CHKERRQ(ierr); 257 ierr = DMInterpolationEvaluate(interpolator, dm, lu, fieldVals);CHKERRQ(ierr); 258 for (p = 0; p < size; ++p) { 259 if (p == rank) { 260 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank);CHKERRQ(ierr); 261 ierr = VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 262 } 263 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 264 } 265 ierr = VecGetArrayRead(interpolator->coords, &vcoords);CHKERRQ(ierr); 266 ierr = VecGetArrayRead(fieldVals, &ivals);CHKERRQ(ierr); 267 for (p = 0; p < interpolator->n; ++p) { 268 for (c = 0; c < Nc; ++c) { 269 #if defined(PETSC_USE_COMPLEX) 270 PetscReal vcoordsReal[3]; 271 PetscInt i; 272 273 for (i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]); 274 #else 275 const PetscReal *vcoordsReal = &vcoords[p*spaceDim]; 276 #endif 277 (*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL); 278 if (PetscAbsScalar(ivals[p*Nc+c] - vals[c]) > PETSC_SQRT_MACHINE_EPSILON) 279 SETERRQ(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); 280 } 281 } 282 ierr = VecRestoreArrayRead(interpolator->coords, &vcoords);CHKERRQ(ierr); 283 ierr = VecRestoreArrayRead(fieldVals, &ivals);CHKERRQ(ierr); 284 /* Cleanup */ 285 ierr = PetscFree(pcoords);CHKERRQ(ierr); 286 ierr = PetscFree2(funcs, vals);CHKERRQ(ierr); 287 ierr = VecDestroy(&fieldVals);CHKERRQ(ierr); 288 ierr = DMRestoreLocalVector(dm, &lu);CHKERRQ(ierr); 289 ierr = DMInterpolationDestroy(&interpolator);CHKERRQ(ierr); 290 ierr = DMDestroy(&dm);CHKERRQ(ierr); 291 ierr = PetscFinalize(); 292 return ierr; 293 } 294 295 /*TEST 296 297 testset: 298 requires: ctetgen 299 args: -dm_plex_dim 3 -petscspace_degree 1 300 301 test: 302 suffix: 0 303 test: 304 suffix: 1 305 args: -dm_refine 2 306 test: 307 suffix: 2 308 nsize: 2 309 args: -petscpartitioner_type simple 310 test: 311 suffix: 3 312 nsize: 2 313 args: -dm_refine 2 -petscpartitioner_type simple 314 test: 315 suffix: 4 316 nsize: 5 317 args: -petscpartitioner_type simple 318 test: 319 suffix: 5 320 nsize: 5 321 args: -dm_refine 2 -petscpartitioner_type simple 322 test: 323 suffix: 6 324 args: -point_type grid 325 test: 326 suffix: 7 327 args: -dm_refine 2 -point_type grid 328 test: 329 suffix: 8 330 nsize: 2 331 args: -petscpartitioner_type simple -point_type grid 332 test: 333 suffix: 9 334 args: -point_type grid_replicated 335 test: 336 suffix: 10 337 nsize: 2 338 args: -petscpartitioner_type simple -point_type grid_replicated 339 test: 340 suffix: 11 341 nsize: 2 342 args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated 343 344 TEST*/ 345