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