1 static char help[] = "Tests for point location\n\n"; 2 3 #include <petscsf.h> 4 #include <petscdmplex.h> 5 6 typedef struct { 7 PetscBool centroids; 8 PetscBool custom; 9 } AppCtx; 10 11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 12 { 13 PetscFunctionBeginUser; 14 options->centroids = PETSC_TRUE; 15 options->custom = PETSC_FALSE; 16 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 } 23 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 } 33 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; 42 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; 55 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 } 73 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; 84 85 PetscFunctionBeginUser; 86 if (!user->custom) PetscFunctionReturn(PETSC_SUCCESS); 87 PetscCall(DMGetCoordinateDim(dm, &cdim)); 88 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)); 100 101 PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF)); 102 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; 113 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)); 125 126 PetscCall(PetscSFDestroy(&cellSF)); 127 PetscCall(VecDestroy(&points)); 128 PetscFunctionReturn(PETSC_SUCCESS); 129 } 130 131 int main(int argc, char **argv) 132 { 133 DM dm; 134 AppCtx user; 135 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 } 146 147 /*TEST 148 149 testset: 150 args: -dm_plex_dim 1 -dm_plex_box_faces 10 151 152 test: 153 suffix: seg 154 155 test: 156 suffix: seg_hash 157 args: -dm_refine 2 -dm_plex_hash_location 158 159 testset: 160 args: -dm_plex_box_faces 5,5 161 162 test: 163 suffix: tri 164 requires: triangle 165 166 test: 167 suffix: tri_hash 168 requires: triangle 169 args: -dm_refine 2 -dm_plex_hash_location 170 171 test: 172 suffix: quad 173 args: -dm_plex_simplex 0 174 175 test: 176 suffix: quad_hash 177 args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location 178 179 testset: 180 args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 181 182 test: 183 suffix: tet 184 requires: ctetgen 185 186 test: 187 suffix: tet_hash 188 requires: ctetgen 189 args: -dm_refine 1 -dm_plex_hash_location 190 191 test: 192 suffix: hex 193 args: -dm_plex_simplex 0 194 195 test: 196 suffix: hex_hash 197 args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location 198 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 203 204 test: 205 suffix: quad_overlap 206 args: -dm_plex_hash_location {{0 1}} 207 208 TEST*/ 209