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 PetscFunctionBegin; 103 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 104 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 105 ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr); 106 ierr = PetscCalloc1(spaceDim,&ind);CHKERRQ(ierr); 107 ierr = PetscCalloc1(spaceDim,&h);CHKERRQ(ierr); 108 h[0] = 1.0/(N-1); h[1] = 1.0/(N-1); h[2] = 1.0/(N-1); 109 ierr = DMDACreate(PetscObjectComm((PetscObject) dm), &da);CHKERRQ(ierr); 110 ierr = DMSetDimension(da, dim);CHKERRQ(ierr); 111 ierr = DMDASetSizes(da, N, N, N);CHKERRQ(ierr); 112 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 113 ierr = DMDASetStencilWidth(da, 1);CHKERRQ(ierr); 114 ierr = DMSetUp(da);CHKERRQ(ierr); 115 ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr); 116 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 117 *Np = info.xm * info.ym * info.zm; 118 ierr = PetscCalloc1(*Np * spaceDim, pcoords);CHKERRQ(ierr); 119 for (k = info.zs; k < info.zs + info.zm; ++k) { 120 ind[2] = k; 121 for (j = info.ys; j < info.ys + info.ym; ++j) { 122 ind[1] = j; 123 for (i = info.xs; i < info.xs + info.xm; ++i, ++n) { 124 ind[0] = i; 125 126 for (d = 0; d < spaceDim; ++d) (*pcoords)[n*spaceDim+d] = ind[d]*h[d]; 127 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D (", rank, n);CHKERRQ(ierr); 128 for (d = 0; d < spaceDim; ++d) { 129 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n*spaceDim+d]);CHKERRQ(ierr); 130 if (d < spaceDim-1) {ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");CHKERRQ(ierr);} 131 } 132 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");CHKERRQ(ierr); 133 } 134 } 135 } 136 ierr = DMDestroy(&da);CHKERRQ(ierr); 137 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr); 138 ierr = PetscFree(ind);CHKERRQ(ierr); 139 ierr = PetscFree(h);CHKERRQ(ierr); 140 *pointsAllProcs = PETSC_FALSE; 141 PetscFunctionReturn(0); 142 } 143 144 static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx) 145 { 146 PetscInt N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d; 147 PetscReal *h; 148 PetscMPIInt rank; 149 PetscErrorCode ierr; 150 151 PetscFunctionBeginUser; 152 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 153 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 154 ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr); 155 ierr = PetscCalloc1(spaceDim,&ind);CHKERRQ(ierr); 156 ierr = PetscCalloc1(spaceDim,&h);CHKERRQ(ierr); 157 h[0] = 1.0/(N-1); h[1] = 1.0/(N-1); h[2] = 1.0/(N-1); 158 *Np = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1); 159 ierr = PetscCalloc1(*Np * spaceDim, pcoords);CHKERRQ(ierr); 160 for (k = 0; k < N; ++k) { 161 ind[2] = k; 162 for (j = 0; j < N; ++j) { 163 ind[1] = j; 164 for (i = 0; i < N; ++i, ++n) { 165 ind[0] = i; 166 167 for (d = 0; d < spaceDim; ++d) (*pcoords)[n*spaceDim+d] = ind[d]*h[d]; 168 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D (", rank, n);CHKERRQ(ierr); 169 for (d = 0; d < spaceDim; ++d) { 170 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n*spaceDim+d]);CHKERRQ(ierr); 171 if (d < spaceDim-1) {ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");CHKERRQ(ierr);} 172 } 173 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");CHKERRQ(ierr); 174 } 175 } 176 } 177 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr); 178 *pointsAllProcs = PETSC_TRUE; 179 ierr = PetscFree(ind);CHKERRQ(ierr); 180 ierr = PetscFree(h);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx) 185 { 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 *pointsAllProcs = PETSC_FALSE; 190 switch (ctx->pointType) { 191 case CENTROID: ierr = CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx);CHKERRQ(ierr);break; 192 case GRID: ierr = CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx);CHKERRQ(ierr);break; 193 case GRID_REPLICATED: ierr = CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx);CHKERRQ(ierr);break; 194 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int) ctx->pointType); 195 } 196 PetscFunctionReturn(0); 197 } 198 199 int main(int argc, char **argv) 200 { 201 AppCtx ctx; 202 PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 203 DM dm; 204 PetscFE fe; 205 DMInterpolationInfo interpolator; 206 Vec lu, fieldVals; 207 PetscScalar *vals; 208 const PetscScalar *ivals, *vcoords; 209 PetscReal *pcoords; 210 PetscBool simplex, pointsAllProcs=PETSC_TRUE; 211 PetscInt dim, spaceDim, c, Np, p; 212 PetscMPIInt rank, size; 213 PetscViewer selfviewer; 214 PetscErrorCode ierr; 215 216 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 217 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 218 ierr = CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);CHKERRQ(ierr); 219 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 220 ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr); 221 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 222 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 223 /* Create points */ 224 ierr = CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx);CHKERRQ(ierr); 225 /* Create interpolator */ 226 ierr = DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator);CHKERRQ(ierr); 227 ierr = DMInterpolationSetDim(interpolator, spaceDim);CHKERRQ(ierr); 228 ierr = DMInterpolationAddPoints(interpolator, Np, pcoords);CHKERRQ(ierr); 229 ierr = DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE);CHKERRQ(ierr); 230 /* Check locations */ 231 for (c = 0; c < interpolator->n; ++c) { 232 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D is in Cell %D\n", rank, c, interpolator->cells[c]);CHKERRQ(ierr); 233 } 234 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr); 235 ierr = VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 236 /* Setup Discretization */ 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, 1, vals, NULL); 278 if (PetscAbsScalar(ivals[p*Nc+c] - vals[c]) > PETSC_SQRT_MACHINE_EPSILON) 279 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); 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: -dm_distribute -petscpartitioner_type simple 310 test: 311 suffix: 3 312 nsize: 2 313 args: -dm_refine 2 -dm_distribute -petscpartitioner_type simple 314 test: 315 suffix: 4 316 nsize: 5 317 args: -dm_distribute -petscpartitioner_type simple 318 test: 319 suffix: 5 320 nsize: 5 321 args: -dm_refine 2 -dm_distribute -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: -dm_distribute -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: -dm_distribute -petscpartitioner_type simple -point_type grid_replicated 339 test: 340 suffix: 11 341 nsize: 2 342 args: -dm_refine 2 -dm_distribute -petscpartitioner_type simple -point_type grid_replicated 343 344 TEST*/ 345