1 static char help[] = "Tests for point location\n\n"; 2 3 #include <petscsf.h> 4 #include <petscdmplex.h> 5 6 /* To inspect the location process, use 7 8 -dm_plex_print_locate 5 -info :dm 9 */ 10 11 typedef struct { 12 PetscBool centroids; 13 PetscBool custom; 14 } AppCtx; 15 16 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 17 { 18 PetscFunctionBeginUser; 19 options->centroids = PETSC_TRUE; 20 options->custom = PETSC_FALSE; 21 22 PetscOptionsBegin(comm, "", "Point Location Options", "DMPLEX"); 23 PetscCall(PetscOptionsBool("-centroids", "Locate cell centroids", "ex17.c", options->centroids, &options->centroids, NULL)); 24 PetscCall(PetscOptionsBool("-custom", "Locate user-defined points", "ex17.c", options->custom, &options->custom, NULL)); 25 PetscOptionsEnd(); 26 PetscFunctionReturn(PETSC_SUCCESS); 27 } 28 29 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) 30 { 31 PetscFunctionBeginUser; 32 PetscCall(DMCreate(comm, dm)); 33 PetscCall(DMSetType(*dm, DMPLEX)); 34 PetscCall(DMSetFromOptions(*dm)); 35 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 36 PetscFunctionReturn(PETSC_SUCCESS); 37 } 38 39 static PetscErrorCode TestCentroidLocation(DM dm, AppCtx *user) 40 { 41 Vec points; 42 PetscSF cellSF = NULL; 43 const PetscSFNode *cells; 44 PetscScalar *a; 45 PetscInt cdim, n; 46 PetscInt cStart, cEnd, c; 47 48 PetscFunctionBeginUser; 49 if (!user->centroids) PetscFunctionReturn(PETSC_SUCCESS); 50 PetscCall(DMGetCoordinateDim(dm, &cdim)); 51 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 52 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 53 /* Locate all centroids */ 54 PetscCall(VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart) * cdim, &points)); 55 PetscCall(VecSetBlockSize(points, cdim)); 56 PetscCall(VecGetArray(points, &a)); 57 for (c = cStart; c < cEnd; ++c) { 58 PetscReal centroid[3]; 59 PetscInt off = (c - cStart) * cdim, d; 60 61 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 62 for (d = 0; d < cdim; ++d) a[off + d] = centroid[d]; 63 } 64 PetscCall(VecRestoreArray(points, &a)); 65 PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF)); 66 PetscCall(VecDestroy(&points)); 67 PetscCall(PetscSFGetGraph(cellSF, NULL, &n, NULL, &cells)); 68 if (n != (cEnd - cStart)) { 69 for (c = 0; c < n; ++c) { 70 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)); 71 } 72 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Located %" PetscInt_FMT " points instead of %" PetscInt_FMT, n, cEnd - cStart); 73 } 74 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); 75 PetscCall(PetscSFDestroy(&cellSF)); 76 PetscFunctionReturn(PETSC_SUCCESS); 77 } 78 79 static PetscErrorCode TestCustomLocation(DM dm, AppCtx *user) 80 { 81 PetscSF cellSF = NULL; 82 const PetscSFNode *cells; 83 const PetscInt *found; 84 Vec points; 85 PetscScalar coords[2] = {0.5, 0.5}; 86 PetscInt cdim, Np = 1, Nfd; 87 PetscMPIInt rank; 88 MPI_Comm comm; 89 90 PetscFunctionBeginUser; 91 if (!user->custom) PetscFunctionReturn(PETSC_SUCCESS); 92 PetscCall(DMGetCoordinateDim(dm, &cdim)); 93 94 // Locate serially on each process 95 PetscCall(VecCreate(PETSC_COMM_SELF, &points)); 96 PetscCall(VecSetBlockSize(points, cdim)); 97 PetscCall(VecSetSizes(points, Np * cdim, PETSC_DETERMINE)); 98 PetscCall(VecSetFromOptions(points)); 99 for (PetscInt p = 0; p < Np; ++p) { 100 const PetscInt idx[2] = {p * cdim, p * cdim + 1}; 101 PetscCall(VecSetValues(points, cdim, idx, coords, INSERT_VALUES)); 102 } 103 PetscCall(VecAssemblyBegin(points)); 104 PetscCall(VecAssemblyEnd(points)); 105 106 PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF)); 107 108 PetscCall(PetscSFGetGraph(cellSF, NULL, &Nfd, &found, &cells)); 109 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 110 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 111 PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found %" PetscInt_FMT " particles\n", rank, Nfd)); 112 for (PetscInt p = 0; p < Nfd; ++p) { 113 const PetscInt point = found ? found[p] : p; 114 const PetscScalar *array; 115 PetscScalar *ccoords = NULL; 116 PetscInt numCoords; 117 PetscBool isDG; 118 119 // Since the v comm is SELF, rank is always 0 120 PetscCall(PetscSynchronizedPrintf(comm, " point %" PetscInt_FMT " cell %" PetscInt_FMT "\n", point, cells[p].index)); 121 PetscCall(DMPlexGetCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords)); 122 for (PetscInt c = 0; c < numCoords / cdim; ++c) { 123 PetscCall(PetscSynchronizedPrintf(comm, " ")); 124 for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscSynchronizedPrintf(comm, " %g", (double)PetscRealPart(ccoords[c * cdim + d]))); 125 PetscCall(PetscSynchronizedPrintf(comm, "\n")); 126 } 127 PetscCall(DMPlexRestoreCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords)); 128 } 129 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 130 131 PetscCall(PetscSFDestroy(&cellSF)); 132 PetscCall(VecDestroy(&points)); 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 136 int main(int argc, char **argv) 137 { 138 DM dm; 139 AppCtx user; 140 141 PetscFunctionBeginUser; 142 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 143 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 144 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 145 PetscCall(TestCentroidLocation(dm, &user)); 146 PetscCall(TestCustomLocation(dm, &user)); 147 PetscCall(DMDestroy(&dm)); 148 PetscCall(PetscFinalize()); 149 return 0; 150 } 151 152 /*TEST 153 154 testset: 155 args: -dm_plex_dim 1 -dm_plex_box_faces 10 156 157 test: 158 suffix: seg 159 160 test: 161 suffix: seg_hash 162 args: -dm_refine 2 -dm_plex_hash_location 163 164 testset: 165 args: -dm_plex_box_faces 5,5 166 167 test: 168 suffix: tri 169 requires: triangle 170 171 test: 172 suffix: tri_hash 173 requires: triangle 174 args: -dm_refine 2 -dm_plex_hash_location 175 176 test: 177 suffix: quad 178 args: -dm_plex_simplex 0 179 180 test: 181 suffix: quad_order_2 182 args: -dm_plex_simplex 0 -dm_coord_petscspace_degree 2 183 184 test: 185 suffix: quad_hash 186 args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location 187 188 testset: 189 args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 190 191 test: 192 suffix: tet 193 requires: ctetgen 194 195 test: 196 suffix: tet_hash 197 requires: ctetgen 198 args: -dm_refine 1 -dm_plex_hash_location 199 200 test: 201 suffix: hex 202 args: -dm_plex_simplex 0 203 204 test: 205 suffix: hex_hash 206 args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location 207 208 test: 209 suffix: hex_order_2 210 args: -dm_plex_simplex 0 -dm_refine 1 -dm_coord_petscspace_degree 2 211 nsize: 2 212 213 test: 214 suffix: hex_order_3 215 args: -dm_plex_simplex 0 -dm_coord_petscspace_degree 3 216 217 testset: 218 args: -centroids 0 -custom \ 219 -dm_plex_simplex 0 -dm_plex_box_faces 21,21 -dm_distribute_overlap 4 -petscpartitioner_type simple 220 nsize: 2 221 222 test: 223 suffix: quad_overlap 224 args: -dm_plex_hash_location {{0 1}} 225 226 # Test location on a Monge Manifold 227 testset: 228 args: -dm_refine 3 -dm_coord_space 0 \ 229 -dm_plex_option_phases proj_ -cdm_proj_dm_plex_coordinate_dim 3 -proj_dm_coord_space \ 230 -proj_dm_coord_remap -proj_dm_coord_map sinusoid -proj_dm_coord_map_params 0.1,1.,1. 231 232 test: 233 requires: triangle 234 suffix: tri_monge 235 236 test: 237 suffix: quad_monge 238 args: -dm_plex_simplex 0 239 240 TEST*/ 241