1c4762a1bSJed Brown static char help[] = "Tests for point location\n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown #include <petscsf.h> 4c4762a1bSJed Brown #include <petscdmplex.h> 5c4762a1bSJed Brown 6*30602db0SMatthew G. Knepley static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) 7c4762a1bSJed Brown { 8c4762a1bSJed Brown PetscErrorCode ierr; 9c4762a1bSJed Brown 10c4762a1bSJed Brown PetscFunctionBeginUser; 11*30602db0SMatthew G. Knepley ierr = DMCreate(comm, dm);CHKERRQ(ierr); 12*30602db0SMatthew G. Knepley ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 13*30602db0SMatthew G. Knepley ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 14c4762a1bSJed Brown ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 15c4762a1bSJed Brown PetscFunctionReturn(0); 16c4762a1bSJed Brown } 17c4762a1bSJed Brown 18*30602db0SMatthew G. Knepley static PetscErrorCode TestLocation(DM dm) 19c4762a1bSJed Brown { 20*30602db0SMatthew G. Knepley PetscInt dim; 21c4762a1bSJed Brown PetscInt cStart, cEnd, c; 22c4762a1bSJed Brown PetscErrorCode ierr; 23c4762a1bSJed Brown 24c4762a1bSJed Brown PetscFunctionBeginUser; 25*30602db0SMatthew G. Knepley ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 26c4762a1bSJed Brown ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 27c4762a1bSJed Brown /* Locate all centroids */ 28c4762a1bSJed Brown for (c = cStart; c < cEnd; ++c) { 29c4762a1bSJed Brown Vec v; 30c4762a1bSJed Brown PetscSF cellSF = NULL; 31c4762a1bSJed Brown const PetscSFNode *cells; 32c4762a1bSJed Brown PetscScalar *a; 33c4762a1bSJed Brown PetscReal centroid[3]; 34c4762a1bSJed Brown PetscInt n, d; 35c4762a1bSJed Brown 36c4762a1bSJed Brown ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 37c4762a1bSJed Brown ierr = VecCreateSeq(PETSC_COMM_SELF, dim, &v);CHKERRQ(ierr); 38c4762a1bSJed Brown ierr = VecSetBlockSize(v, dim);CHKERRQ(ierr); 39c4762a1bSJed Brown ierr = VecGetArray(v, &a);CHKERRQ(ierr); 40c4762a1bSJed Brown for (d = 0; d < dim; ++d) a[d] = centroid[d]; 41c4762a1bSJed Brown ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 42c4762a1bSJed Brown ierr = DMLocatePoints(dm, v, DM_POINTLOCATION_NONE, &cellSF);CHKERRQ(ierr); 43c4762a1bSJed Brown ierr = VecDestroy(&v);CHKERRQ(ierr); 44c4762a1bSJed Brown ierr = PetscSFGetGraph(cellSF,NULL,&n,NULL,&cells);CHKERRQ(ierr); 45c4762a1bSJed Brown if (n != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Found %d cells instead %d", n, 1); 46c4762a1bSJed Brown 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); 47c4762a1bSJed Brown ierr = PetscSFDestroy(&cellSF);CHKERRQ(ierr); 48c4762a1bSJed Brown } 49c4762a1bSJed Brown PetscFunctionReturn(0); 50c4762a1bSJed Brown } 51c4762a1bSJed Brown 52c4762a1bSJed Brown int main(int argc, char **argv) 53c4762a1bSJed Brown { 54c4762a1bSJed Brown DM dm; 55c4762a1bSJed Brown PetscErrorCode ierr; 56c4762a1bSJed Brown 57c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 58*30602db0SMatthew G. Knepley ierr = CreateMesh(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr); 59*30602db0SMatthew G. Knepley ierr = TestLocation(dm);CHKERRQ(ierr); 60c4762a1bSJed Brown ierr = DMDestroy(&dm);CHKERRQ(ierr); 61c4762a1bSJed Brown ierr = PetscFinalize(); 62c4762a1bSJed Brown return ierr; 63c4762a1bSJed Brown } 64c4762a1bSJed Brown 65c4762a1bSJed Brown /*TEST 66c4762a1bSJed Brown 67c4762a1bSJed Brown test: 68c4762a1bSJed Brown suffix: 0 69c4762a1bSJed Brown requires: triangle 70*30602db0SMatthew G. Knepley args: -dm_coord_space 0 -dm_view ascii::ascii_info_detail 71c4762a1bSJed Brown 72c4762a1bSJed Brown TEST*/ 73