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