1 static char help[] = "Tests for point location\n\n"; 2 3 #include <petscsf.h> 4 #include <petscdmplex.h> 5 6 typedef struct { 7 /* Domain and mesh definition */ 8 PetscInt dim; /* The topological mesh dimension */ 9 PetscBool cellSimplex; /* Use simplices or hexes */ 10 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ 11 PetscBool testPartition; /* Use a fixed partitioning for testing */ 12 PetscInt testNum; /* Labels the different test partitions */ 13 } AppCtx; 14 15 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 16 { 17 PetscErrorCode ierr; 18 19 PetscFunctionBeginUser; 20 options->dim = 2; 21 options->cellSimplex = PETSC_TRUE; 22 options->filename[0] = '\0'; 23 options->testPartition = PETSC_TRUE; 24 options->testNum = 0; 25 26 ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 27 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex13.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 28 ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex13.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsString("-filename", "The mesh file", "ex13.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex13.c", options->testPartition, &options->testPartition, NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsBoundedInt("-test_num", "The test partition number", "ex13.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 32 ierr = PetscOptionsEnd(); 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 37 { 38 DM dmDist = NULL; 39 PetscInt dim = user->dim; 40 PetscBool cellSimplex = user->cellSimplex; 41 const char *filename = user->filename; 42 size_t len; 43 PetscErrorCode ierr; 44 45 PetscFunctionBeginUser; 46 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 47 if (len) {ierr = DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);CHKERRQ(ierr);} 48 else {ierr = DMPlexCreateBoxMesh(comm, dim, cellSimplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);} 49 if (user->testPartition) { 50 PetscPartitioner part; 51 PetscInt *sizes = NULL; 52 PetscInt *points = NULL; 53 PetscMPIInt rank, size; 54 55 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 56 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 57 if (!rank) { 58 if (dim == 2 && cellSimplex && size == 2) { 59 switch (user->testNum) { 60 case 0: { 61 PetscInt triSizes_p2[2] = {4, 4}; 62 PetscInt triPoints_p2[8] = {3, 5, 6, 7, 0, 1, 2, 4}; 63 64 ierr = PetscMalloc2(2, &sizes, 8, &points);CHKERRQ(ierr); 65 ierr = PetscArraycpy(sizes, triSizes_p2, 2);CHKERRQ(ierr); 66 ierr = PetscArraycpy(points, triPoints_p2, 8);CHKERRQ(ierr); 67 break;} 68 case 1: { 69 PetscInt triSizes_p2[2] = {6, 2}; 70 PetscInt triPoints_p2[8] = {1, 2, 3, 4, 6, 7, 0, 5}; 71 72 ierr = PetscMalloc2(2, &sizes, 8, &points);CHKERRQ(ierr); 73 ierr = PetscArraycpy(sizes, triSizes_p2, 2);CHKERRQ(ierr); 74 ierr = PetscArraycpy(points, triPoints_p2, 8);CHKERRQ(ierr); 75 break;} 76 default: 77 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 78 } 79 } else if (dim == 2 && cellSimplex && size == 3) { 80 PetscInt triSizes_p3[3] = {3, 3, 2}; 81 PetscInt triPoints_p3[8] = {1, 2, 4, 3, 6, 7, 0, 5}; 82 83 ierr = PetscMalloc2(3, &sizes, 8, &points);CHKERRQ(ierr); 84 ierr = PetscArraycpy(sizes, triSizes_p3, 3);CHKERRQ(ierr); 85 ierr = PetscArraycpy(points, triPoints_p3, 8);CHKERRQ(ierr); 86 } else if (dim == 2 && !cellSimplex && size == 2) { 87 PetscInt quadSizes_p2[2] = {2, 2}; 88 PetscInt quadPoints_p2[4] = {2, 3, 0, 1}; 89 90 ierr = PetscMalloc2(2, &sizes, 4, &points);CHKERRQ(ierr); 91 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 92 ierr = PetscArraycpy(points, quadPoints_p2, 4);CHKERRQ(ierr); 93 } else SETERRQ3(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition dim: %D simplex: %D size: %D", dim, (PetscInt) cellSimplex, (PetscInt) size); 94 } 95 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 96 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr); 97 ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr); 98 ierr = PetscFree2(sizes, points);CHKERRQ(ierr); 99 } 100 ierr = DMPlexDistribute(*dm, 0, NULL, &dmDist);CHKERRQ(ierr); 101 if (dmDist) { 102 ierr = DMDestroy(dm);CHKERRQ(ierr); 103 *dm = dmDist; 104 } 105 ierr = PetscObjectSetName((PetscObject) *dm, cellSimplex ? "Simplicial Mesh" : "Tensor Product Mesh");CHKERRQ(ierr); 106 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode TestLocation(DM dm, AppCtx *user) 111 { 112 PetscInt dim = user->dim; 113 PetscInt cStart, cEnd, c; 114 PetscErrorCode ierr; 115 116 PetscFunctionBeginUser; 117 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 118 /* Locate all centroids */ 119 for (c = cStart; c < cEnd; ++c) { 120 Vec v; 121 PetscSF cellSF = NULL; 122 const PetscSFNode *cells; 123 PetscScalar *a; 124 PetscReal centroid[3]; 125 PetscInt n, d; 126 127 ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 128 ierr = VecCreateSeq(PETSC_COMM_SELF, dim, &v);CHKERRQ(ierr); 129 ierr = VecSetBlockSize(v, dim);CHKERRQ(ierr); 130 ierr = VecGetArray(v, &a);CHKERRQ(ierr); 131 for (d = 0; d < dim; ++d) a[d] = centroid[d]; 132 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 133 ierr = DMLocatePoints(dm, v, DM_POINTLOCATION_NONE, &cellSF);CHKERRQ(ierr); 134 ierr = VecDestroy(&v);CHKERRQ(ierr); 135 ierr = PetscSFGetGraph(cellSF,NULL,&n,NULL,&cells);CHKERRQ(ierr); 136 if (n != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Found %d cells instead %d", n, 1); 137 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); 138 ierr = PetscSFDestroy(&cellSF);CHKERRQ(ierr); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 int main(int argc, char **argv) 144 { 145 DM dm; 146 AppCtx user; /* user-defined work context */ 147 PetscErrorCode ierr; 148 149 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 150 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 151 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 152 ierr = TestLocation(dm, &user);CHKERRQ(ierr); 153 ierr = DMDestroy(&dm);CHKERRQ(ierr); 154 ierr = PetscFinalize(); 155 return ierr; 156 } 157 158 /*TEST 159 160 test: 161 suffix: 0 162 requires: triangle 163 args: -test_partition 0 -dm_view ascii::ascii_info_detail 164 165 TEST*/ 166