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(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 276 for (PetscInt i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]); 277 #else 278 const PetscReal *vcoordsReal = &vcoords[p * spaceDim]; 279 #endif 280 PetscCall((*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL)); 281 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); 282 } 283 } 284 PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords)); 285 PetscCall(VecRestoreArrayRead(fieldVals, &ivals)); 286 /* Cleanup */ 287 PetscCall(PetscFree(pcoords)); 288 PetscCall(PetscFree2(funcs, vals)); 289 PetscCall(VecDestroy(&fieldVals)); 290 PetscCall(DMRestoreLocalVector(dm, &lu)); 291 PetscCall(DMInterpolationDestroy(&interpolator)); 292 PetscCall(DMDestroy(&dm)); 293 PetscCall(PetscFinalize()); 294 return 0; 295 } 296 297 /*TEST 298 299 testset: 300 requires: ctetgen 301 args: -dm_plex_dim 3 -petscspace_degree 1 302 303 test: 304 suffix: 0 305 test: 306 suffix: 1 307 args: -dm_refine 2 308 test: 309 suffix: 2 310 nsize: 2 311 args: -petscpartitioner_type simple 312 test: 313 suffix: 3 314 nsize: 2 315 args: -dm_refine 2 -petscpartitioner_type simple 316 test: 317 suffix: 4 318 nsize: 5 319 args: -petscpartitioner_type simple 320 test: 321 suffix: 5 322 nsize: 5 323 args: -dm_refine 2 -petscpartitioner_type simple 324 test: 325 suffix: 6 326 args: -point_type grid 327 test: 328 suffix: 7 329 args: -dm_refine 2 -point_type grid 330 test: 331 suffix: 8 332 nsize: 2 333 args: -petscpartitioner_type simple -point_type grid 334 test: 335 suffix: 9 336 args: -point_type grid_replicated 337 test: 338 suffix: 10 339 nsize: 2 340 args: -petscpartitioner_type simple -point_type grid_replicated 341 test: 342 suffix: 11 343 nsize: 2 344 args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated 345 346 TEST*/ 347