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