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