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