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