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 output_file: output/empty.out 157 158 test: 159 suffix: seg 160 161 test: 162 suffix: seg_hash 163 args: -dm_refine 2 -dm_plex_hash_location 164 165 testset: 166 args: -dm_plex_box_faces 5,5 167 output_file: output/empty.out 168 169 test: 170 suffix: tri 171 requires: triangle 172 173 test: 174 suffix: tri_hash 175 requires: triangle 176 args: -dm_refine 2 -dm_plex_hash_location 177 178 test: 179 suffix: quad 180 args: -dm_plex_simplex 0 181 182 test: 183 suffix: quad_order_2 184 args: -dm_plex_simplex 0 -dm_coord_petscspace_degree 2 185 186 test: 187 suffix: quad_hash 188 args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location 189 190 testset: 191 args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 192 output_file: output/empty.out 193 194 test: 195 suffix: tet 196 requires: ctetgen 197 198 test: 199 suffix: tet_hash 200 requires: ctetgen 201 args: -dm_refine 1 -dm_plex_hash_location 202 203 test: 204 suffix: hex 205 args: -dm_plex_simplex 0 206 207 test: 208 suffix: hex_hash 209 args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location 210 211 test: 212 suffix: hex_order_2 213 args: -dm_plex_simplex 0 -dm_refine 1 -dm_coord_petscspace_degree 2 214 nsize: 2 215 216 test: 217 suffix: hex_order_3 218 args: -dm_plex_simplex 0 -dm_coord_petscspace_degree 3 219 220 testset: 221 args: -centroids 0 -custom \ 222 -dm_plex_simplex 0 -dm_plex_box_faces 21,21 -dm_distribute_overlap 4 -petscpartitioner_type simple 223 nsize: 2 224 225 test: 226 suffix: quad_overlap 227 args: -dm_plex_hash_location {{0 1}} 228 229 # Test location on a Monge Manifold 230 testset: 231 args: -dm_refine 3 -dm_coord_space 0 \ 232 -dm_plex_option_phases proj_ -cdm_proj_dm_plex_coordinate_dim 3 -proj_dm_coord_space \ 233 -proj_dm_coord_remap -proj_dm_coord_map sinusoid -proj_dm_coord_map_params 0.1,1.,1. 234 output_file: output/empty.out 235 236 test: 237 requires: triangle 238 suffix: tri_monge 239 240 test: 241 suffix: quad_monge 242 args: -dm_plex_simplex 0 243 244 TEST*/ 245