Actual source code: ex17.c

  1: static char help[] = "Tests for point location\n\n";

  3: #include <petscsf.h>
  4: #include <petscdmplex.h>

  6: typedef struct {
  7:   PetscBool centroids;
  8:   PetscBool custom;
  9: } AppCtx;

 11: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 12: {
 13:   PetscFunctionBeginUser;
 14:   options->centroids = PETSC_TRUE;
 15:   options->custom    = PETSC_FALSE;

 17:   PetscOptionsBegin(comm, "", "Point Location Options", "DMPLEX");
 18:   PetscCall(PetscOptionsBool("-centroids", "Locate cell centroids", "ex17.c", options->centroids, &options->centroids, NULL));
 19:   PetscCall(PetscOptionsBool("-custom", "Locate user-defined points", "ex17.c", options->custom, &options->custom, NULL));
 20:   PetscOptionsEnd();
 21:   PetscFunctionReturn(PETSC_SUCCESS);
 22: }

 24: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
 25: {
 26:   PetscFunctionBeginUser;
 27:   PetscCall(DMCreate(comm, dm));
 28:   PetscCall(DMSetType(*dm, DMPLEX));
 29:   PetscCall(DMSetFromOptions(*dm));
 30:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: static PetscErrorCode TestCentroidLocation(DM dm, AppCtx *user)
 35: {
 36:   Vec                points;
 37:   PetscSF            cellSF = NULL;
 38:   const PetscSFNode *cells;
 39:   PetscScalar       *a;
 40:   PetscInt           cdim, n;
 41:   PetscInt           cStart, cEnd, c;

 43:   PetscFunctionBeginUser;
 44:   if (!user->centroids) PetscFunctionReturn(PETSC_SUCCESS);
 45:   PetscCall(DMGetCoordinateDim(dm, &cdim));
 46:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
 47:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 48:   /* Locate all centroids */
 49:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart) * cdim, &points));
 50:   PetscCall(VecSetBlockSize(points, cdim));
 51:   PetscCall(VecGetArray(points, &a));
 52:   for (c = cStart; c < cEnd; ++c) {
 53:     PetscReal centroid[3];
 54:     PetscInt  off = (c - cStart) * cdim, d;

 56:     PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
 57:     for (d = 0; d < cdim; ++d) a[off + d] = centroid[d];
 58:   }
 59:   PetscCall(VecRestoreArray(points, &a));
 60:   PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF));
 61:   PetscCall(VecDestroy(&points));
 62:   PetscCall(PetscSFGetGraph(cellSF, NULL, &n, NULL, &cells));
 63:   if (n != (cEnd - cStart)) {
 64:     for (c = 0; c < n; ++c) {
 65:       if (cells[c].index != c + cStart) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Could not locate centroid of cell %" PetscInt_FMT ", error %" PetscInt_FMT "\n", c + cStart, cells[c].index));
 66:     }
 67:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Located %" PetscInt_FMT " points instead of %" PetscInt_FMT, n, cEnd - cStart);
 68:   }
 69:   for (c = cStart; c < cEnd; ++c) PetscCheck(cells[c - cStart].index == c, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate centroid of cell %" PetscInt_FMT ", instead found %" PetscInt_FMT, c, cells[c - cStart].index);
 70:   PetscCall(PetscSFDestroy(&cellSF));
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

 74: static PetscErrorCode TestCustomLocation(DM dm, AppCtx *user)
 75: {
 76:   PetscSF            cellSF = NULL;
 77:   const PetscSFNode *cells;
 78:   const PetscInt    *found;
 79:   Vec                points;
 80:   PetscScalar        coords[2] = {0.5, 0.5};
 81:   PetscInt           cdim, Np = 1, Nfd;
 82:   PetscMPIInt        rank;
 83:   MPI_Comm           comm;

 85:   PetscFunctionBeginUser;
 86:   if (!user->custom) PetscFunctionReturn(PETSC_SUCCESS);
 87:   PetscCall(DMGetCoordinateDim(dm, &cdim));

 89:   // Locate serially on each process
 90:   PetscCall(VecCreate(PETSC_COMM_SELF, &points));
 91:   PetscCall(VecSetBlockSize(points, cdim));
 92:   PetscCall(VecSetSizes(points, Np * cdim, PETSC_DETERMINE));
 93:   PetscCall(VecSetFromOptions(points));
 94:   for (PetscInt p = 0; p < Np; ++p) {
 95:     const PetscInt idx[2] = {p * cdim, p * cdim + 1};
 96:     PetscCall(VecSetValues(points, cdim, idx, coords, INSERT_VALUES));
 97:   }
 98:   PetscCall(VecAssemblyBegin(points));
 99:   PetscCall(VecAssemblyEnd(points));

101:   PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF));

103:   PetscCall(PetscSFGetGraph(cellSF, NULL, &Nfd, &found, &cells));
104:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
105:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
106:   PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found %" PetscInt_FMT " particles\n", rank, Nfd));
107:   for (PetscInt p = 0; p < Nfd; ++p) {
108:     const PetscInt     point = found ? found[p] : p;
109:     const PetscScalar *array;
110:     PetscScalar       *ccoords = NULL;
111:     PetscInt           numCoords;
112:     PetscBool          isDG;

114:     // Since the v comm is SELF, rank is always 0
115:     PetscCall(PetscSynchronizedPrintf(comm, "  point %" PetscInt_FMT " cell %" PetscInt_FMT "\n", point, cells[p].index));
116:     PetscCall(DMPlexGetCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords));
117:     for (PetscInt c = 0; c < numCoords / cdim; ++c) {
118:       PetscCall(PetscSynchronizedPrintf(comm, "  "));
119:       for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscSynchronizedPrintf(comm, " %g", (double)PetscRealPart(ccoords[c * cdim + d])));
120:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
121:     }
122:     PetscCall(DMPlexRestoreCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords));
123:   }
124:   PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));

126:   PetscCall(PetscSFDestroy(&cellSF));
127:   PetscCall(VecDestroy(&points));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: int main(int argc, char **argv)
132: {
133:   DM     dm;
134:   AppCtx user;

136:   PetscFunctionBeginUser;
137:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
138:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
139:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
140:   PetscCall(TestCentroidLocation(dm, &user));
141:   PetscCall(TestCustomLocation(dm, &user));
142:   PetscCall(DMDestroy(&dm));
143:   PetscCall(PetscFinalize());
144:   return 0;
145: }

147: /*TEST

149:   testset:
150:     args: -dm_plex_dim 1 -dm_plex_box_faces 10

152:     test:
153:       suffix: seg

155:     test:
156:       suffix: seg_hash
157:       args: -dm_refine 2 -dm_plex_hash_location

159:   testset:
160:     args: -dm_plex_box_faces 5,5

162:     test:
163:       suffix: tri
164:       requires: triangle

166:     test:
167:       suffix: tri_hash
168:       requires: triangle
169:       args: -dm_refine 2 -dm_plex_hash_location

171:     test:
172:       suffix: quad
173:       args: -dm_plex_simplex 0

175:     test:
176:       suffix: quad_hash
177:       args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location

179:   testset:
180:     args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3

182:     test:
183:       suffix: tet
184:       requires: ctetgen

186:     test:
187:       suffix: tet_hash
188:       requires: ctetgen
189:       args: -dm_refine 1 -dm_plex_hash_location

191:     test:
192:       suffix: hex
193:       args: -dm_plex_simplex 0

195:     test:
196:       suffix: hex_hash
197:       args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location

199:   testset:
200:     args: -centroids 0 -custom \
201:           -dm_plex_simplex 0 -dm_plex_box_faces 21,21 -dm_distribute_overlap 4 -petscpartitioner_type simple
202:     nsize: 2

204:     test:
205:       suffix: quad_overlap
206:       args: -dm_plex_hash_location {{0 1}}

208: TEST*/