1 static char help[] = "Tests for point location\n\n"; 2 3 #include <petscsf.h> 4 #include <petscdmplex.h> 5 6 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) { 7 PetscFunctionBeginUser; 8 PetscCall(DMCreate(comm, dm)); 9 PetscCall(DMSetType(*dm, DMPLEX)); 10 PetscCall(DMSetFromOptions(*dm)); 11 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 12 PetscFunctionReturn(0); 13 } 14 15 static PetscErrorCode TestLocation(DM dm) { 16 Vec points; 17 PetscSF cellSF = NULL; 18 const PetscSFNode *cells; 19 PetscScalar *a; 20 PetscInt cdim, n; 21 PetscInt cStart, cEnd, c; 22 23 PetscFunctionBeginUser; 24 PetscCall(DMGetCoordinateDim(dm, &cdim)); 25 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 26 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 27 /* Locate all centroids */ 28 PetscCall(VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart) * cdim, &points)); 29 PetscCall(VecSetBlockSize(points, cdim)); 30 PetscCall(VecGetArray(points, &a)); 31 for (c = cStart; c < cEnd; ++c) { 32 PetscReal centroid[3]; 33 PetscInt off = (c - cStart) * cdim, d; 34 35 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 36 for (d = 0; d < cdim; ++d) a[off + d] = centroid[d]; 37 } 38 PetscCall(VecRestoreArray(points, &a)); 39 PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF)); 40 PetscCall(VecDestroy(&points)); 41 PetscCall(PetscSFGetGraph(cellSF, NULL, &n, NULL, &cells)); 42 if (n != (cEnd - cStart)) { 43 for (c = 0; c < n; ++c) { 44 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)); 45 } 46 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Located %" PetscInt_FMT " points instead of %" PetscInt_FMT, n, cEnd - cStart); 47 } 48 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); 49 PetscCall(PetscSFDestroy(&cellSF)); 50 PetscFunctionReturn(0); 51 } 52 53 int main(int argc, char **argv) { 54 DM dm; 55 56 PetscFunctionBeginUser; 57 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 58 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 59 PetscCall(TestLocation(dm)); 60 PetscCall(DMDestroy(&dm)); 61 PetscCall(PetscFinalize()); 62 return 0; 63 } 64 65 /*TEST 66 67 testset: 68 args: -dm_plex_dim 1 -dm_plex_box_faces 10 69 70 test: 71 suffix: seg 72 73 test: 74 suffix: seg_hash 75 args: -dm_refine 2 -dm_plex_hash_location 76 77 testset: 78 args: -dm_plex_box_faces 5,5 79 80 test: 81 suffix: tri 82 requires: triangle 83 84 test: 85 suffix: tri_hash 86 requires: triangle 87 args: -dm_refine 2 -dm_plex_hash_location 88 89 test: 90 suffix: quad 91 args: -dm_plex_simplex 0 92 93 test: 94 suffix: quad_hash 95 args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location 96 97 testset: 98 args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 99 100 test: 101 suffix: tet 102 requires: ctetgen 103 104 test: 105 suffix: tet_hash 106 requires: ctetgen 107 args: -dm_refine 1 -dm_plex_hash_location 108 109 test: 110 suffix: hex 111 args: -dm_plex_simplex 0 112 113 test: 114 suffix: hex_hash 115 args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location 116 117 TEST*/ 118