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 PetscErrorCode ierr; 9 10 PetscFunctionBeginUser; 11 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 12 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 13 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 14 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 static PetscErrorCode TestLocation(DM dm) 19 { 20 PetscInt dim; 21 PetscInt cStart, cEnd, c; 22 PetscErrorCode ierr; 23 24 PetscFunctionBeginUser; 25 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 26 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 27 /* Locate all centroids */ 28 for (c = cStart; c < cEnd; ++c) { 29 Vec v; 30 PetscSF cellSF = NULL; 31 const PetscSFNode *cells; 32 PetscScalar *a; 33 PetscReal centroid[3]; 34 PetscInt n, d; 35 36 ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 37 ierr = VecCreateSeq(PETSC_COMM_SELF, dim, &v);CHKERRQ(ierr); 38 ierr = VecSetBlockSize(v, dim);CHKERRQ(ierr); 39 ierr = VecGetArray(v, &a);CHKERRQ(ierr); 40 for (d = 0; d < dim; ++d) a[d] = centroid[d]; 41 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 42 ierr = DMLocatePoints(dm, v, DM_POINTLOCATION_NONE, &cellSF);CHKERRQ(ierr); 43 ierr = VecDestroy(&v);CHKERRQ(ierr); 44 ierr = PetscSFGetGraph(cellSF,NULL,&n,NULL,&cells);CHKERRQ(ierr); 45 if (n != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D: Found %d cells instead of 1", c, n); 46 if (cells[0].index != c) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate centroid of cell %D, instead found %D", c, cells[0].index); 47 ierr = PetscSFDestroy(&cellSF);CHKERRQ(ierr); 48 } 49 PetscFunctionReturn(0); 50 } 51 52 int main(int argc, char **argv) 53 { 54 DM dm; 55 PetscErrorCode ierr; 56 57 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 58 ierr = CreateMesh(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr); 59 ierr = TestLocation(dm);CHKERRQ(ierr); 60 ierr = DMDestroy(&dm);CHKERRQ(ierr); 61 ierr = PetscFinalize(); 62 return ierr; 63 } 64 65 /*TEST 66 67 test: 68 suffix: seg 69 args: -dm_plex_dim 1 -dm_plex_box_faces 10 70 71 test: 72 suffix: tri 73 requires: triangle 74 args: -dm_coord_space 0 -dm_plex_box_faces 5,5 75 76 test: 77 suffix: quad 78 args: -dm_plex_simplex 0 -dm_plex_box_faces 5,5 79 80 test: 81 suffix: tet 82 requires: ctetgen 83 args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 84 85 test: 86 suffix: hex 87 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 88 89 TEST*/ 90